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ABSTRACT 


In  this  thesis  we  study  several  computer  implementa- 
tions of  the  thinning  algorithm,  a  new  method  for  generating 
non-homogeneous  Poisson  processes.   The  method,  developed 
by  Professor  P.A.W.  Lewis,  Naval  Postgraduate  School, 
Monterey,  California,  and  G.S.  Shedler,  IBM  Research 
Laboratory,  San  Jose,,  California,  is  valid  for  Poisson 
processes  with  any  given  intensity  function.   The  basic 
thinning  algorithm  is  modified  to  exploit  several  refine- 
ments which  reduce  computer  execution  time  by  approximately 
one-third.   The  basic  and  modified  thinning  programs  are 
compared  with  a  previous  algorithm  of  Lewis  and  Shedler, 
the  Poisson  decomposition  and  gap-statistics  algorithm, 

which  is  easily  implemented  for  Poisson  processes  with 

2 

intensity  tunctions  of  the  form  exp (aQ+a, t+a^t  ).   The 

thinning  programs  are  competitive  in  both  execution  time 
and  computer  memory  requirements.   One  program  implementa- 
tion generates  the  events  in  a  Poisson  process  one  at  a 
time;  another  program  im^plements  the  algorithmic  refinements 
which  improve  efficiency. 
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I.   INTRODUCTION 

The  Poisson  process  is  a  widely  known  and  studied 
stochastic  process.   It  is  frequently  used  to  describe 
random  arrivals  at  some  type  of  service  facility  such  as 
a  service  station  fuel  pump  or  a  bank  teller's  windov/. 
In  its  most  common  form,  the  "rate"  of  these  arrivals  is 
considered  to  be  constant  over  time.   This  is  the  homogeneous 
Poisson  process  which  has  the  familiar  property  that  times 
between  arrivals  (or  events)  are  exponentially  distributed 
with  mean  equal  to  the  inverse  of  the  rate. 

The  assumption  of  a  constant  rate,  or  homogeneity,  is 
at  best  tenuous  when  applied  to  real  world  data.   For  exam- 
ple, the  rate  of  arrivals  at  a  traffic  light  typically 
varies  from  very  high  during  rush  hours  to  very  low  in  the 
early  morning.   In  addition  to  this  cyclic  time-of-day 
effect,  arrival  rates  may  exhibit  longer  term  increases  or 
decreases.   Further,  these  effects  may  be  superimposed  upon 
shorter  term  effects  to  produce  a  more  complex  rate  which 
varies  with  time.   These  processes  for  which  arrival  rates 
vary  with  time  may  often  be  represented  by  a  non-homogeneous 
Poisson  process,  that  is,  a  Poisson  process  with  a  time 
dependent  rate  of  arrival. 

The  generic  term  of  Poisson  process  includes  then  both 
homogeneous  and  non-homogeneous  Poisson  processes.   LEWIS 
and  SHEDLER  [Ref.  1]  define  the  Poisson  process  generally  in 
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terms  of  a  monotone  non-decreasing  right-continous  function 
A(t)  which  is  bounded  in  any  finite  interval.   Then  the 
number  of  points,  N(t'',t'),  in  any  finite  interval  has  a 
Poisson  distribution  with  parameter  y(t'',t')  =  A(t')  -  A(t'') 
Thus,  for  example  in  (0,t']/  with  t'  >_  0,    P{N(t'',t')  =  n} 
E  P{N^,=n}  =  pQ  e    /nl,  where  Uq  =  y(0,t')  =  A(t')  -  A(0). 

The  right  derivative  a (t)  of  A(t)  will  be  assumed  to 
exist  and  is  called  the  rate  function  or  intensity  function 
of  the  process.   A(t)  is  called  the  integrated  rate  function 
and  has  the  interpretation  that  for  t  ^  0,  A(t)  -  A(0)  =  E[N  ] 
For  the  homogeneous  Poisson  process  X (t)  is  a  constant, 
e.g.  A,  and  thus  the  integrated  rate  function  is  simply  the 
product  of  A  and  t,  i.e.  the  expected  value  of  N  =   N(0,t). 

While  simulation  of  homogeneous  Poisson  processes  is 
relatively  straightforward,  the  non-homogeneous  Poisson 
process  is  more  problematical.   Times  between  events  are 
not  exponential  in  the  general  case  and  simulation  has 
typically  been  tailored  to  specific  classes  of  intensity 
functions.   LEWIS  a-nd  SHEDLER  [Ref  •  1]  list  three  general 
methods  for  simulating  non-homogeneous  Poisson  processes 
and  one  method  for  a  special  rate  function.   The  general 
methods  include  the  time  scale  transformation  method  and 
the  conditioning  and  order-statistics  method.   The  special 
method  is  the  gap-statistics  method,  a  method  which  is 
particular  to  the  degree-one  exponential  polynomial  intensity 
fucntion,  i.e.  those  of  the  form  A (t)  =  explb^  +  b^t) . 
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Implementation  of  the  general  methods  on  a  computer  may 
pose  special  problems.   Often  the  inverse  of  the  integrated 
rate  function  is  not  explicit  and  must  be  computed  numerically, 
Other  problems  in  implementation  generally  result  in  lower 
efficiency,  as  measured  by  execution  time  or  computer  storage 
requirements  or  both. 

One  class  of  intensity  functions  which  is  of  general 
interest  is  the  degree-two  exponential  polynomial  family. 

That  is,  those  with  intensity  function  of  the  form 

2 

A(t)  =  exp(a„  +  a,t  +  a„t  ).   This  family  of  functions  has 

the  property  of  being  positive  for  all  values  of  t,  a 
necessary  condition  for  an  intensity  function.   Additionally, 
by  varying  the  magnitude  and  sign  of  the  coefficients,  the 
exponential  polynomial  of  degree  two  can  be  made  to  be  mono- 
tone increasing  or  decreasing  over  time,  as  well  as  increasing 
and  then  decreasing,  or  vice  versa.   Use  of  this  type  of 
intensity  function  also  leads  to  statistical  procedures 
which  are  relatively  simple. 

LEWIS  and  SHEDLER  [Ref.  2]  proposed  a  new  method  of 
generating  the  non-homogeneous  Poisson  process  with  degree- 
two  exponential  polynomial  intensity  function.   It  involves 
decomposition  of  the  degree-two  exponential  polynomial 
intensity  function,  X(t),  into  two  functions,  a  degree-one 
exponential  polynomial  function,  A  (t) ,  and  a  difference 
function,  a_(-c)  =  X  (t)  -  X  (t)  .   This  procedure  allows  the 

D  J-i 

points  in  the  degree-one  exponential  polynomial  event  stream 
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to  be  generated  using  the  gap-statistics  method,  which  is 
highly  efficient  when  implemented  on  a  computer.   The 
remaining  points  with  intensity  function  ^^(t)  are  generated 
by  other  methods  and  then  merged  with  the  other  events. 

PATROW  [Ref.  3]  implemented  two  algorithms,  the  time 
scale  transformation  algorithm  and  the  Poisson  decomposi- 
tion and  gap-statistics  technique,  and  compared  them  for 
computational  speed  and  computer  memory  requirements.   His 
results  indicated  that  the  Poisson  decomposition  and  gap- 
statistics  technique  was  from  two  to  seven  times  faster  than 
the  time  scale  transformation  algorithm,  although  the  former 
required  about  thirty  percent  more  computer  memory. 

PATROW s  work  [Ref.  3]  is  also  an  excellent  self- 
contained  reference  on  Poisson  processes,  bringing  many 
references  together  under  one  cover. 

A  recent  result  of  LEWIS  and  SHEDLER  [Ref.  1]  develops 
a  new  method  for  generation  of  points  in  a  non-homogeneous 
Poisson  process.   This  method,  called  "thinning",  is  similar 
to  the  general  conditioning-acceptance-rejection  method 
but  has  subtle  differences  which  are  computationally  signi- 
ficant.  The  thinning  method  is  straightforward  in  both  an 
analytical  and  a  computational  sense,  and  is  valid  for  any 
type  of  intensity  function.   The  thinning  theorem  is 
presented  in  Section  II. 

This  thesis  is,  in  a  sense,  a  sequel  to  PATROW s  work 
[Ref.  3].   Its  purpose  is  to  implement  the  thinning  algorithm 
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in  computer  program  form  and  to  compare  it  to  the  Poisson 
decomposition  and  gap-statistics  algorithm  implemented  by 
PATROW  [Ref.  3].   The  latter  implementation  was  designed 
for  a  specific  subset  of  intensity  functions,  degree-two 
exponential  polynomials.   Since  the  Poisson  decomposition 
and  gap-statistics  method  outperformed  a  general  case 
algorithm  (time  scale  transformation)  by  a  significant 
margin,  comparing  the  thinning  method  to  the  Poisson 
decomposition  and  gap-statistics  method  should  give  a 
reasonable  indication  of  the  thinning  algorithm's  performance 
in  generating  non-homogeneous  Poisson  processes  with  other 
than  degree-two  exponential  polynomial  intensity  functions. 
Section  III  lists  the  two  algorithms  considered, as  well 
as  a  special  application  of  the  thinning  process  which  will 
be  of  interest  to  those  involved  in  event-step  simulation. 
Section  IV  describes  the  methodology  used  in  comparing  the 
algorithms  while  Section  V  deals  with  aspects  of  the  thinning 
procedure  which  may  be  exploited  to  enhance  its  overall 
effectiveness  in  a  variety  of  situations.   Finally,  Section 
VI  presents  the  results  and  conclusions  of  the  comparisons 
of  the  algorithms.   Appendices  A  and  B  contain  secondary 
results  and  computer  program  listings  following  the 
appendices . 


II.   THE  THINNING  THEOREM 

The  underlying  concept  of  the  thinning  method  involves 

* 

the  use  of  a  "bounding"  Poisson  process,  {N   :  t  >_  0 }  , 

* 
where  N   is  the  number  of  points  in  the  bounding  process  in 

the  interval  (0,t].   This  process  may  be  either  homogeneous 

or  non-homogeneous  Poisson,  but  should  be  one  which  is  easy 

to  simulate  on  a  computer.   It  is  called  bounding  because 

its  intensity  function,  denoted  X    (t) ,  bounds  the  intensity 

function  A(t),  of  the  nonhomogeneous  Poisson  process  which 

is  to  be  simulated  over  the  fixed  interval  (0,t'].   That 

it  "k 

is,  X    (t)  ^  A(t)  for  all  t  in  (0,t'].   Points  at  T., 
i  =  1,  .../  N  , ,  are  generated  for  the  bounding  process 
over  the  interval  (0,t']-   These  points  are  then  deleted, 
or  "thinned",  with  independent  probabilities  equal  to 

1  -  {X{T.)/X    (T.)).   Thus  the  probability  that  a  point  of 

* 

the  bounding  process,  T.,  is  a  point  of  the  process  being 

generated  is  equal  to  the  ratio  of  the  intensity  functions 

*      ie  -k 

evaluated  at  that  point,  i.e.  a(T.)/X  (T. ) . 
More  formally: 

Theorem  1.   Consider  the  one-dimensional  non-homogeneous 

*  * 

Poisson  process  {N   :  t  >_  0 }  with  rate  function  a  (t)  . 

The  number  of  events,  N  , ,  in  the  fixed  interval  (0,t'] 
has  a  Poisson  distribution  with  parameter  y  (0,t')  =  y 
=  A*(t')  -  A*(0) . 
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*    *    *  * 

Let  T,,  T^/  T^,  .../  T  *  ,  be  the  times  of  the  events 

of  the  process  in  the  interval  (O/t']- 

Suppose  that  for  0  <  t  <  t ' ,  A (t)  <  A  (t) .   For 


i  =  1 


,  2,  .../  N.  ,  delete  the  event  at  T.  with  independent 


probability  1  -  A(T.)/A  (T.). 

Then  the  remaining  times  form  a  non-homogeneous  Poisson 
process  with  rate  function  A(t)  in  the  interval  (0,t']. 
Proof: 

We  assume  that  A(t)  is  continuous  and  use  the  definition 
of  the  Poisson  process  based  on  incremental  probabilities. 
Thus  we  need  to  show  that  the  occurrence  of  an  event  in 
(t/t+dt]  is  independent  of  the  number  or  times  of  occurrence 
of  events  before  t,  and  that 


^^\+dt  -  ^t  "  °^  "  ^  "  '^^^^^^  "^  o(dt). 


P^^N^^^^  -  N^  =  1}  =  A(t)dt  +  o(dt), 
t+dt     t 


and 


P^^t-fdt  -  \  >  i>  =  °(dt; 


Now  we  have  that 
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Pino  event  from  {n   :  t  >_  0 }  in  (t,t+dt]} 

* 

=  P{no  event  from  (N   :  t  ^  0}  in  (t,t+dt]}+  P{ event 

* 

from  {N^  :  t  >  0}  in  (t,t+dt]  and  it  is  "thinned"} 

=  1  -  X    (t)dt  +  [A  (t)dt] • [1  -  A(t)/X*(t) ]  +  o (dt) 


=  1  -  A  (t)dt  +  A  (t)dt  -  A*(t)  -A  (t)/A* (t)  -dt  +  o (dt! 


=  1  -  A (t)dt  +  o(dt) 


Similarly: 


P{one  event  from  {N   :  t  ^  0 }  in  (t,t+dtj} 

* 

=  P{event  from  (N   :  t  ^  0}  in  (t,t+dt]  which  is  not  "thinned"} 


=  A  (t)dt'A(t)/A  (t)  +  o(dt) 
=  A  (t)dt  +  o(t) 

Also  it  follows  directly  that 

P{more  than  one  event  in  (t,t+dt]}   =   o(dt) 

Moreover,  an  event  in  (t,t+dt]  is  independent  of  what 
happens  before  t  because: 
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1.  {N   :  t  >_  0}  is  a  Poisson  process  and  therefore  has 
independent  increments ,  and 

2.  The  thinning  uniform  random  variate  is  independent 

of  other  thinning  variates,  and  is  independent  of  the 

* 

Poisson  process  {N   :  t  >_  0 } . 

Q.E.D. 

Figure  1  shows  a  graphical  representation  of  a  particu- 
lar case  of  bounding  and  object  intensity  functions. 
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III.   ALGORITHMS  CONSIDERED 

A.   POISSON  DECOMPOSITION  AND  GAP-STATISTICS  ALGORITHM 
1.   Usage 

This  algorithm  is  the  one  found  by  PATROW  [Ref.  3] 
to  be  most  efficient  in  simulation  of  the  degree-two  exponen- 
tial polynomial  class  of  intensity  functions  and  its  imple- 
mentation by  PATROW  was  confined  to  that  group.   Basically, 

the  approach  is  to  decompose  the  intensity  function  (which 

2 
IS  of  the  form  A(t)  =  exp(a^  +  a,t  +  a^t  ))  into  a  degree- 
one  exponential  polynomial  function,  A  (t)  =  exp(b„  +  b-t), 

L  U      X 

and  a  difference  function,  Ar)(t)  =  A  (t)  -  A  (t)  .   The  points 
or  events  in  the  process  with  the  degree-one  exponential 
polynomial  function,  A  (t) ,  are  generated  over  the  interval 

Lt 

(0,t']  utilizing  the  computationally  fast  gap-statistics 
algorithm.   The  points  in  the  process  with  intensity  func- 
tion A_^(t)  are  generated  using  conditioning-acceptance- 
rejection  techniques.   The  two  event  streams  are  then 
superposed  to  produce  the  event  stream  for  the  non-homogeneous 
Poisson  process  with  the  intensity  function  A (t) . 

In  the  case  where  A(t)  has  an  internal  maximum 
or  minimum  in  the  interval  (0,t']/  the  interval  is  par- 
titioned and  treated  as  two  separate  intervals  for  simula- 
tion with  the  event  streams  being  merged  in  the  final  step. 

Efficiency  is  optimized  by  maximizing  the  area 
under  the  degree-one  exponential  polynomial  intensity 
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function,  A  (t) ,  subject,  of  course,  to  the  constraint 
A^ (t)  <  A(t)  for  0  <  t  <  t'.   This  maximizes  the  use  of 
the  faster  gap-statistics  algorithm  and  minimizes  the  use 
of  the  conditioning-acceptance-rejection  algorithm  which 
is  slow  relative  to  the  gap-statistics  algorithm. 

PATROW  [Ref.  3]  deals  extensively  with  the  details 
of  this  algorithm  and  it  is  consequently  presented  here 
only  in  outline  form. 

2 .   Algorithm  Statement 

a.  Categorize  the  intensity  function,  A (t) ,  into 

one  of  six  cases  by  examination  of  the  coefficients  a, 

2 
and  a„  in  A (t)  =  exp(aQ  +  a,t  +  a^t  ).   Examples  of  each 

of  these  cases  are  shown  in  Figures  2  through  Figure  7 

located  in  Section  IV. 

b.  (1)   If  A (t)  is  monotone  increasing  or  monotone 
decreasing  over  the  interval  (Cases  I,  II,  IV  and  V;  see 
Figures  2,3,5  and  6),  decompose  A (t)  into  A  (t) ,  which  is 
degree-one  exponential  polynomial,  and  A  (t)  =  A (t)  -  A  (t) 
Thus  the  decomposed  functions  have  the  forms : 


A^(t)   =   exp(bQ  +  b^t) 


Aj^(t)   =   exp(aQ  +  a^t  +  a^t^)  -  exp(bQ  +  b^t) 


and 


A(t)   =   A^(t)  +  Aj^(t) 
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Choose  bj,  and  b,  so  as  to  maximize  the 
area  under  A  (t)  subject  to  A^ (t)  <  A(t)  for  all  t  in 

Li  J_i        — 

(0,t'] . 

(2)   If  A(t)  is  not  monotone  over  the  interval 

(Cases  III  and  VI;  see  Figures  4  and  7) ,  partition  the 
interval  (0,t']  into  two  disjoint,  contiguous  subintervals , 

(0/b]  and  (b,t']-   Choose  b  as  the  (unique)  point  where 
A (b)  is  a  maximum  (minimum)  of  A (t)  over  (0,t'].   Treat 
each  subinterval  as  in  b.(l),  applying  subsequent  steps  on 
each  subinterval  separately,  and  combining  results  as  the 
final  step. 

c.  Generate  points  in  the  Poisson  process  with 
degree-one  exponential  polynomial  intensity  function, 

A  (t) ,  using  the  gap-statistics  method. 

d.  Generate  and  order  points  in  the  Poisson  process 
with  intensity  function  A  (t)  using  the  conditioning- 
acceptance-rejection  method. 

e.  Merge  (superpose)  the  two  event  streams  from 
Step  3  and  Step  4.   The  merged  stream  is  from  the  non- 
homogeneous  Poisson  process  with  intensity  function  A(t). 

B.   THE  BASIC  THINNING  ALGORITHM 
1.   Usage 

The  thinning  theorem  is  implemented  in  a  straight- 
forward manner.   Typically,  the  bounding  process  used  is 

*  * 

homogeneous  Poisson  with  constant  rate  A  ,  where  A   is  an 
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upper  bound  of  A(t)  over  (0,t'].  In  this  case  efficiency 
is  optimized  if  A  is  the  least  upper  bound  (LUB)  of  A (t) 
over  (0,t']. 

2.   Algorithm  Statement 

a.  Generate  events  in  the  Poisson  process 

*  * 

{N^  :  t  ^  0}  with  rate  function  A  (t)  in  the  fixed  interval 

(0,t'].   If  the  number  of  events  generated,  n  ,  is  such 

* 

that  n   =0,  exit;  there  are  no  events  in  the  process 

{N   :  t  ^  0}. 

b.  Denote  the  (ordered)  events  by  T, ,  T^ ,  ...,  T  *. 

12        n 

Set  i  =  1  and  k  =  0 . 

c.  Generate  U.,  uniformly  distributed  between  0 

and  1.   If  U.  <_  A(T.)/A  (T.),  set  k  equal  to  k+1  and  T,  =  T. 

* 

d.  Set  i  equal  to  i+1.   If  i  <_  n  ,  go  to  c. 

e.  Return  T, ,  T^ ,  ...,  T  ,  where  n  =  k,  and  also  n. 

12         n 

C.   THE  ONE-AT-A-TIME  THINNING  ALGORITHM 
1.   Usage 

In  some  event-step  simulations,  it  is  customary  or 
necessary  to  generate  only  one  event  at  a  time,  rather  than 
an  array  of  events.   The  thinning  algorithm  is  easily 
modified  to  generate  the  next  event  in  the  non-homogeneous 
Poisson  process  with  intensity  function  A (t) .   In  this  case, 
the  algorithm  utilizes  the  time  of  the  last  event,  T-_-i; 
the  right  hand  limit,  t',  of  the  fixed  interval  over  which 

the  process  is  being  simulated;  and  the  bounding  process 

* 

intensity  function,  A  (t) .   All  variates  are  generated 
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one  at  a  time,  thus  no  arrays  are  required  for  storage. 

The  output  is  T. ,  the  time  of  the  next  event,  if  any,  in 

the  interval  (T.  wt']- 

1-1 

The  algorithm  is  stated  here  for  the  case  in  which 

the  bounding  process  is  homogeneous  Poisson,  i.e., 

*       * 
X    (t)  =  A  ,  a  constant  and  an  upper  bound  of  X (t) . 

Specifically: 

T-  is  obtained  by  generating  and  cumulating 

*  *      * 

exponential  (mean  =  1/A  )  random  variates  E.  , ,  E.  „,  ... 

1 ,  J.    1  /  ^ 

for  i  =  1,  2,  ...,  until  for  the  first  time. 


•  *     * 

U.   .<A(T.  ,+E.  ,+...+E.   .)/A 
i/D  -     1-1     1/1  1/3 


A  detailed  algorithmic  statement  of  this  procedure 

follows. 

2 .   Algorithm  Statement 

If  i  =  1,  set  T.  ,  =  0  (i.e.  the  left  end  point 

1-1 

of  the  interval)  ,  otheirwise,  T.  ,  is  known.   Then  for  each 

1-1 

i  =  1,  2,  ...  ,  the  time,  T.,  of  the  event  in  the  non- 
homogeneous  Poisson  process  is  given  by  the  following: 

a.  Set  j  =  1 

* 

b.  Generate  E.  . ,  an  exponential  random  variate 

i/D    . 

with  mean  1/A  .   If  T.  ,  +   V  E-  ,  is  greater  than  t', 

1-1    k=l  "-'^ 
exit;  there  are  no  more  points  in  the  interval  (T._,,t']. 

c.  Generate  U.,  uniformly  distributed  between  0 
and  1.   If 
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J 

k=l 


set 


* 

T.   =   T.  ,  +   )   E.  , 
k=l 

and  exit. 

d.   Otherwise  set  j  =  j+1;  go  to  b. 

Note:   U.  .  and  U.  are  uniformly  distributed 
between  0  and  1. 
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IV.   METHODOLOGY  FOR  ALGORITHM  COMPARISON 

A.   MEASURES  OF  EFFECTIVENESS 

Two  quantifiable  measures  of  effectiveness  were  chosen 
\        as  yardsticks  for  algorithm  comparison.   These  were  compu- 
tational speed  and  computer  memory  requirements.   Som.e  other 
considerations,  such  as  programming  ease  and  robustness, 
are  discussed  in  Section  VI.   It  must  also  be  recalled  that 
the  classes  of  intensity  functions  for  which  the  two 
algorithms  are  usable  are  different.   The  Poisson  decomposi- 
tion and  gap-statistics  algorithm  is  only  easily  implemented 

for  a  restricted  set  of  intensity  functions,  those  of  the 

2 
form  X (t)  =  exp(a„  +  a,t  +  a^t  ),  i.e.  the  degree-two 

exponential  polynomial.   Conversely,  the  thinning  algorithm 

is  valid  for  any  positive  intensity  function.   Thus  a  direct 

comparison  can  only  be  made  in  that  subset  of  intensity 

functions  for  which  both  algorithm  implementations  are 

valid,  the  degree- two  exponential  polynomials. 

PATROW  [Ref.  3]  developed  six  sample  intensity  functions, 
all  special  cases  of  the  degree-two  exponential  polynomial, 
and  these  are  used  herein  as  the  test  cases.   These  are 
described  in  Section  IV. C. 3  below. 

1.   Computational  Speed 

Typically,  computer  time  is  a  costly  commodity 
in  economic  terms.   It  may  also  be  a  significant  factor 
in  determining  more  mundane  considerations  such  as  job 
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priority  and  thus  turn-around  time.  Thus  computer  run 
time  is  a  natural  candidate  as  a  measure  of  effectiveness 
for  comparing  competing  algorithms. 

PATROW  [Ref.  3]  utilized  a  procedure  in  which 
event  streams  from  each  of  six  sample  intensity  functions 
were  replicated  several  times  in  "packages".   The  number 
of  replications  was  large  if  the  expected  number  of  events 
in  the  event  stream  was  small,  and  vice  versa.   Thus  the 
product  of  the  number  of  events  times  the  number  of  replica- 
tions was  kept  on  the  same  order  of  magnitude.   For  sim.plicity, 
the  same  technique  was  used  here  although  results  showed  a 
wide  variation  in  the  run  times  for  the  six  packages. 

Programming  of  the  thinning  algorithm  was  done  so  as 
to  minimize  run  time  while  maintaining  parity  with  the 
Poisson  decomposition  and  gap-statistics  algorithm  wherever 
direct  comparison  could  be  made.   For  example,  shuffled 
random  numbers  are  called  in  both  programs  and  both  are 
dimensioned  to  accommodate  event  streams  of  up  to  5000  events. 

Undoubtedly  further  programming  refinements  exist 
which  might  increase  slightly  the  speed  of  one  or  the  other 
algorithm.   Also,  different  computers  might  have  unique 
features  which  could  be  exploited.   The  overall  purpose  here 
was  to  obtain  a  relative  order  of  magnitude  comparison  and 
it  believed  that  this  objective  was  accomplished  in  every 
meaningful  sense. 
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2.   Computer  Memory  Requirements 

This  is  the  second  obvious  means  of  comparing  two 
algorithms.   Again,  some  core  reduction  could  undoubtedly 
be  made  by  a  sophisticated  programmer.   Most  notably,  core 
requirements  can  be  reduced  substantially  if  only  one  non- 
homogeneous  Poisson  variate  is  generated  at  a  time  (the 
one-at-a-time  algorithm)  but  this  has  the  predictable 
effect  of  increasing  execution  time  considerably  (see 
Tables  IV,  V,  and  VI) . 

B.   MEASUREMENT  CONSIDERATIONS 

Measurement  of  computer  memory  requirements  is  straight- 
foJTward  and  deterministic. 

Measurement  of  computational  speed,  more  specifically 
Central  Processing  Unit  (CPU)  time,  is  quite  another  matter. 
First  of  all,  the  number  of  events  in  each  replication  of 
the  non-homogeneous  Poisson  process  varies  causing  CPU  time 
to  be  a  random  variable.   More  important,  however,  are  the 
effects  of  internal  computer  procedures. 

In  the  first  place,  the  so-called  CPU  time  printed  out 
on  the  normal  IBM-360  output  has  only  a  general  relationship 
to  the  actual  computational  time  required  by  the  CPU.   This 
is  caused  by  the  addition  of  certain  "overhead"  time.   This 
overhead  time  is  a  function  of  the  number  of  other  programs 
in  the  system  as  well  as  such  factors  as  compilation  and 
linkage  times.   Thus  the  same  program,  run  at  two  different 
times  may  differ  in  "CPU  time"  by  a  factor  of  two  or  more. 


Program  execution  times  were  isolated  from  compilation 
and  linkage  times  by  the  use  of  a  system  subroutine, 
GETIME.   This  subroutine  allows  the  user  to  initialize  an 
internal  timer  within  the  program  and  read  cumulative  time 
at  various  points  in  the  program.   Although  this  method  is 
not  exact,  it  does  measure  actual  elapsed  CPU  time  to  within 
a  small  fraction  of  a  second.   This  does  not,  however, 
entirely  alleviate  the  time-of-day  effect  experienced  when 
running  the  same  program  at  different  times.   That  is, 
although  the  elapsed  CPU  time  can  be  measured  accurately, 
the  same  program  will  generally  have  somewhat  different 
execution  times  each  time  it  is  run  [Ref.  4].   Theoretically, 
the  execution  times  would  be  constant  for  stand-alone  runs, 
i.e.  runs  with  no  other  competing  programs  in  the  system. 
This  is  rarely  realized  in  practice. 

These  considerations  lead  to  the  development  and  use  of 
the  side-by-side  setup  described  below.   This  method  appears 
to  be  statistically  sound  as  a  means  of  dealing  with  the 
problems  of  time  measurement.   Due  to  the  differences  in 
execution  times  noted,  the  best  measure  of  effectiveness 
was  determined  to  be  a  ratio  of  execution  times  for  the 
respective  algorithms,  rather  than  absolute  times. 

C.   TEST  SETUP 

1.   Computational  Speed 

The  central  idea  here  was  to  equalize  the  effects 
of  non-essential  processes  on  each  algorithm.   This  was 
accomplished  by  the  following  algorithm: 
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1.  Set  k  =  1. 

2.  Zero  internal  time  clock. 

3.  Call  Algorithm  A.   Replicate  M  times. 

4.  Read  internal  time  clock.   Store  time. 

5.  Zero  internal  time  clock. 

6.  Call  Algorithm  B.   Replicate  M  times. 

7.  Read  internal  time  clock.   Store  time. 

8.  Set  k  =  k  +  1.   If  k  is  greater  than  k    ,  go  to  9 . 

^  max   ^ 

Otherwise  go  to  2. 

9 .  Compute  m.ean  and  variance  of  the  k     execution  times 

^  max 

for  each  algorithm. 

10.   Compute  ratio  of  means. 

This  procedure  was  used  in  all  comparisons.   M, 

the  number  of  replications  per  package,  varied  between  30 

and  100  as  discussed  above.   K    ,  the  number  of  times  each 

max 

package  was  replicated,  was  typically  set  equal  to  thirty. 
2 .   Computer  Memory  Requirements 

To  measure  computer  memory  requirements,  a  small 
main  program  was  written,  calling  the  subroutine  which 
implemented  the  program  being  measured.   Total  program  length 
in  bytes  was  obtained  from  the  standard  computer  output  and 
the  core  alloted  to  the  main  program  was  subtracted  to 
obtain  the  desired  figure.   This  includes  all  library  rou- 
tines and  arrays  for  storage  of  event  times  and  arrays  of 
random  variates. 
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The  core  requirements  are  deterministic  in  that 
they  do  not  change  from  one  run  to  another  but  are  strictly 
a  function  of  the  program  coding. 
3.   Test  Cases  Utilized 

PATROW  [Ref.  3]  developed  six  sample  intensity 
functions  representing  the  possible  variations  in  sign  and 

relative  magnitude  of  the  coefficients  a-,  and  a^  in  the 

2 

exponential  polynomial,  exp(a^  +  a,t  +  a^t  ). 

Since  the  sam.ple  intensity  functions  were  designed 
to  test  different  aspects  of  the  Poisson  decomposition 
and  gap-statistics  algorithm,  they  were  also  used  for  com- 
parison here.   Although  each  algorithm  is  affected  by 
different  considerations,  the  test  cases  do,  coincidentally , 
put  the  thinning  algorithm  through  its  paces. 

For  continuity,  the  six  test  cases,  or  sample 
intensity  functions,  are  presented  below  in  Figures  2  through 
7. 
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V.   EFFICIENCY  AND  PROGRAMMING  CONSIDERATIONS 

A.   GENERAL 

This  section  deals  with  factors  which  affect  the  per- 
formance of  the  thinning  algorithm.   Four  specific  areas 
are  presented  in  which  significant  gains  in  terms  of  com- 
putational speed  may  be  realized.   With  the  exception  of 
Section  V.D  (which  applies  only  to  the  case  of  exponential 
polynomial  intensity  functions)  these  considerations  apply 
to  the  general  class  of  intensity  functions. 

In  general  application,  one  of  the  primary  indicators 
of  efficiency  is  the  relative  size  of  the  area  under  the 
intensity  function  to  that  of  the  bounding  function,  i.e. 
the  ratio,  R,  given  by: 


f         f 
R   =   /   A(t)dt//   X*(t)dt 
0  0 


Since  both  numerator  and  denominator  are  simply  the  respec- 
tive integrated  rate  functions,  A(t),  evaluated  at  either 

end  of  the  interval,  R  is  the  ratio  of  the  expected  number 

* 

of  events  in  the  two  processes,  i.e.  E[N  , ]/E[N  , ]. 

Case  1  of  the  sample  intensity  functions  is  particularly 

illustrative  (see  Figure  2) .   The  intensity  function 

2 
A(t)  =  exp(1.6  +  O.OlSt  +  O.OOOSt  )  is  bounded  on  the 
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interval  (0,100]  by  a  least  upper  bound  (LUB)  of  3294.47. 

If  a  homogeneous  Poisson  process  with  rate  equal  to  the 

* 
LUB  is  used  as  the  bounding  function,  E[N  , ]  =  329,447 

points  will  be  generated  on  the  average.   Of  these,  all 
but  1464  will  be  rejected  on  the  average  (i.e.  E [N  , ] ) . 
The  ratio  of  the  respective  expected  values  is  thus 
1464/329,447  =  0.0044  =  the  ratio  of  the  areas  under  the 
intensity  functions. 

Thus  a  rough  relative  measure  of  the  efficiency  of  the 
thinning  process  in  a  particular  situation  can  be  gained  by 
examining  a  graph  of  the  two  intensity  functions,  even  if 
the  expected  values  are  not  easily  calculated.   This  pro- 
cedure may  also  be  an  indicator  in  deciding  whether  to 
partition  the  interval  and  use  different  bounding  functions 
on  each  subinterval. 

B.   UTILIZATION  OF  ARRAYS  OF  RANDOM  VARIATES 

Computer  generated  random  variates  are  used  both  in 
generating  the  points  of  the  bounding  process  {N   :  t  >_  0 } 
and  in  the  actual  thinning  process  itself.   Since  the 
number  of  variates  required  is  typically  large,  efficient 
generation  becomes  a  programming  consideration  for  medium 
to  large  scale  simulations. 

The  basic  thinning  program  presented  herein  requires 
both  exponential  and  uniform  random  variates .   Both  types 
are  obtained  utilizing  the  random  number  package,  LLRANDOM, 
developed  by  LEARMONTH  and  LEWIS  [Ref.  5].   Shuffled 
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exponential  random  variates  of  mean  1.0  are  generated  using 
the  SEXPON  subroutine  while  shuffled  uniform  (0,1)  random 
variates  are  obtained  from  the  SRAND  subroutine.   Both  of 
these  routines  offer  considerable  "economies  of  scale"  in 
terms  of  time  when  multiple  numbers  or  arrays  of  variates 
are  generated  at  once,  as  opposed  to  one-at-a-time  genera- 
tion.  Using  the  test  setup  of  Section  IV. C,  average  times 
to  generate  varying  quantities  of  random  numbers  were 
determined.   Table  I  reveals  the  relative  savings  realized 
by  calling  large  arrays  of  random  numbers.   Thus  considerable 
time  can  be  saved  by  generating  all  required  random  num.bers 
from  one  subroutine  call.   Programming  difficulties  involve 
deciding  how  many  variates  to  generate.   The  general  goal 
is  to  generate  as  many  as  needed  while  keeping  the  unused 
excess  to  a  minimiim.   The  balance  used  was  to  generate  the 
expected  number  required  plus  an  excess  of  four  standard 
deviations.   For  example,  in  the  generation  of  the  bounding 
process,  the  expected  number  of  points,  E[N  , ]  is  A  -t' 
and  the  variance  is  the  same.   Thus  the  number  of  exponen- 

r-  *  ... 

tials  called  was  y  +  4/y  where  y  =  A  -t'.   Provision  is 
made  for  the  unlikely  (1  in  40,000)  case  that  more  are 
required. 

For  specific  applications  this  procedure  could  be 
improved  slightly.   For  example,  if  the  expected  number  of 
points,  E [N  , ] ,  is  small,  e.g.  100,  then  the  expected 
excess  (four  standard  deviations)  comprises  forty  percent 
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Type  of 
Variate 


Number 
Called 


Total  Time   Mean  Time  Per 
(ysec)      Variate  (ysec) 


Exponential 

1 

784 

784 

Exponential 

10 

1293 

129 

Exponential 

100 

7343 

73 

Exponential 

1000 

68046 

68 

Uniform 

1 

1213 

1213 

Uniform 

10 

1381 

138 

Uniform 

100 

3276 

33 

Uniform 

1000 

21544 

22 

Sample  Size 

=  200  (each 
Table  I 

groupin 

Generation  Times  for  Arrays  of  Shuffled  Random  Variates 

From  LLRANDOM 


of  the  total  whereas  for  large  E [N  , ] ,  e.g.  4000,  the 
expected  "waste"  is  only  about  six  percent.   In  the  former 
case,  reducing  the  "padding"  to  one  or  two  standard  devia- 
tions would,  on  the  average,  increase  efficiency  slightly 
although  the  probability  of  a  second  subroutine  call  for 
more  random  variates  would  increase. 

As  an  example,  if  we  were  to  call  100  exponential 
variates  one  at  a  time,  the  total  time  is  78,400  usee, 
compared  to  7,343  ysec  for  100  exponential  variates  called 


in  an  array.   For  100  +  4/100  =  140  variates,  the  time 
is  still  small  compared  to  78,400  usee. 
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C.   UTILIZATION  OF  INTENSITY  FUNCTION  LOWER  BOUND 

One  of  the  most  time  consuming  repetitive  operations 
is  the  computation  of  the  intensity  function  value,  A(t), 
during  the  thinning  process.   In  the  case  of  the  exponential 
polynomial  intensity  function,  this  involves  one  power,  two 
multiplications,  two  additions,  and  one  exponentiation  for 
each  point  generated  in  the  bounding  process.   Since  points 
are  accepted  for  the  nor>-.homogeneous  Poisson  process  when- 
ever the  uniform  (0,1)  thinning  variate  is  less  than  the 
ratio  A(t)/A  (t) ,  considerable  time  savings  result  if  the 
intensity  function  has  a  positive  lower  bound,  say  A_, 
since  points  are  always  accepted  when  the  uniform  (0,1) 
variate  is  less  than  the  ratio  X/\    .   In  the  general  case, 
this  ratio  must  be  calculated  only  once.   The  expected  num- 
ber of  intensity  function  computations  which  are  alleviated 

*     * 

by  the  use  of  the  lower  bound  is  given  by  (_a/A  )  E  [N,  ,  ]  where 

A_  is  a  lower  bound  of  the  intensity  function;  \   is  an  upper 

bound  of  the  intensity  function  (both  bounds  are  over  the 

* 
interval  (0,t'])  and  E[N  , ]  is  the  average  number  of  points 

to  be  thinned,  i.e.  the  average  number  of  events  in  the 

bounding  process. 

It  is  clear  that  the  closer  A_  is  to  being  the  greatest 

* 
lower  bound  (GLB)  and  the  closer  X      is  to  being  the  LUB, 

the  more  efficient  the  program. 

If  the  intensity  function  is  strictly  non-decreasing 

a  further  (and  potentially  great)  improvement  is  realized 

by  initially  setting  A  equal  to  A,  and  then  setting  it 
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subsequently  equal  to  the  last  value  of  the  intensity 
function,  A (t) .   This  results  in  a  monotone  increasing 
lower  bound  and  thus  a  decreasing  probability  of  evaluating 
the  intensity  function. 

Test  cases  II  through  VI  were  run  side  by  side  with  and 
without  the  use  of  a  lower  bound  for  the  intensity  function. 
On  the  average,  the  program  which  did  not  utilize  a  lower 
bound  required  twenty  percent  more  time  than  the  program 
using  a  lower  bound.   Please  see  Appendix  B  for  case-by- 
case  comparison. 

D.   UTILIZATION  OF  EXPONENTIAL  VARIATES  FOR  THINNING  OF 
EXPONENTIAL  POLYNOMIAL  INTENSITY  FUNCTIONS 

The  time  requirements  for  evaluating  A (t)  were  discussed 

in  Section  V.C  above.   In  the  case  of  exponential  polynomial 

2 
intensity  functions,  e.g.  A (t)  =  exp(a^  +  a,t  +  a„t  ),  the 

major  contributor  to  computation  time  is  the  exponentiation 

operation.   Exponentiation  can  be  avoided  by  utilizing  the 

following  relationship: 


* 

U.  <_   A(t)/A      if  and  only  if 


■*  *  9 

E.   =   -In  U.  >_  InA   -  InA  (t)  =  InA   -  (a^  +  a^t  +  a^t    ), 


where : 


U.  is  a  uniform  (0,1)  random  variate, 

E.  is  an  exponential  random  variate  with  m.ean  one 
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Thus  the  thinning  test  to  accept  points  from  the 
bounding  process  becomes: 


*  2  * 

If  E^  ^  InX   -  (a^  +  a^t  +  a2t  ) ,    for  t  =  T. ,  accept 

* 

T.  as  a  point  in  the  non- homogeneous  Poisson  process 

with  rate  A(t);  othei>;ise,  reject  T.  (i.e.  thin  it). 


The  key  to  this  relationship  lies  in  the  fact  that  if 
U  is  distributed  uniform  (0,1),  then  -In  U  is  distributed 
as  a  unit  exponential  variate,  i.e.  an  exponential  variate 
with  mean  one.   This  is  shown  by  the  following: 

Let  U  be  uniform  (0,1). 

Then  P{U  ^  x}   E   p{ln  U  £  In  x}   e   P{-ln  U  >  -In  x} 

but  P{U  £  x}   =   X,  thus  let  y   =   -In  x, 

then  P{-ln  U  _>  y}   =   exp(-y). 

Thus  -In  U  is  distributed  as  a  unit  exponential  variate. 

Although  more  time  is  required  to  generate  the  exponential 
random  variates  for  thinning  than  the  uniforms,  the  alle- 
viation of  the  exponentiation  operation  more  than  compen- 
sates for  the  additional  generation  time.   This  is  because 
SEXPON,  the  portion  of  LLRANDOM  which  generates  exponentials, 
generates  exponential  variates  by  the  Marsaglia  "rectangle- 
wedge-triangle"  method,  which  is  faster  than  taking  logarithms 

Since  exponential  random  variates  are  used  in  the  genera- 
tion of  the  bounding  homogeneous  Poisson  process,  an  addi- 
tional time  savings  can  be  realized  by  using  the  variates 
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which  are  left  over  (i.e.  not  used)  from  generating  the 
bounding  process  (these  are  generated  in  arrays) . 

For  the  test  cases  considered,  use  of  exponentials 
for  thinning  resulted  in  an  average  time  savings  of  ten 
percent.   Please  see  Appendix  B  for  case-by-case  results. 

E.   RECYCLING  OF  THINNING  VARIATES 

As  mentioned  above,  a  uniform  or  exponential  random 
variate  is  required  for  each  point  to  be  "thinned" .   Each 
of  these  variates  requires  a  significant  amount  of  time  for 
generation.   Obviously  a  time  savings  would  be  realized  if 
fewer  variates  were  required. 

1.   Recycling  of  Uniform  Variates 

Assume  U.  is  uniform  (0,1)  but  that  its  value  is 
1 

unknown.   Assume  then  that  further  information  becomes 
available  that  U.  is  less  than  a  (0  <  a  <  1),  but  its  value 
is  still  unknown.   Then  U.  is  uniformly  distributed  over 
the  interval  (0,a).   If  U.  ,  is  then  computed  by  "scaling 
up"  U.,  i.e.  dividing  U.  by  a,  then  U.  ,  is  uniform  (0,1). 
Similarly,  if  U.  is  uniform  (0,1)  and  subsequent  infojrmation 
places  it  somewhere  above  a,  then  U.  ,  =  (U.  -  a)/(l  -  a) 
is  uniform  (0,1).   Thus  by  conditioning  on  whether  the 
variate  is  greater  than  or  less  than  a  given  value,  a  new 
variate  can  be  computed  with  the  desired  properties. 
Moreover,  this  variate  is  independent  of  its  predecessor. 

In  the  thinning  algorithm,  each  point  is  tested 

* 

using  a  uniform  (0,1)  variate.   Specifically,  if  U^  <_  A(T^)/A 
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the  point  T.  is  accepted  as  a  point  in  the  non-homogeneous 
Poisson  process.   Since  the  ratio  X(T.)/A   is  between  zero 
and  one,  and  the  only  test  is  whether  U.  is  less  than  or 

■k  -k 

greater  than  the  ratio  A(T^)/A  ,  the  next  uniform  (0,1) 
variate,  ^j_+i'    ^^^  ^^  generated  using  the  rules  above. 
The  algorithm  is: 

1.  Let  U.  be  uniform  (0,1).   If  U.  is  less  than 

1  1 

a  =  A(T*)/A*,  let  U^^^  =  U^/ ( A (T*) /A*l ;  exit. 

2.  Otherwise  let  U^^^  =  [U^  -  (A(T^)/A  ] / [1- ( A (T . ) /A  )] 

U.,,  is  uniform  (0,1). 
1+1 

In  theory,  only  one  uniform  random  variate  is 
required  for  the  entire  thinning  process  I   In  computational 
practice,  however,  care  must  be  exercised  because  of  the 
finite  capacity  of  the  computer  to  represent  numbers.   After 
ten  to  twenty  divisions  the  scaled  uniform  number  will 
consist  only  of  low-order  bits  of  the  random  number  and 
these  are  usually  not  uniformly  distributed. 

If  the  intensity  function  has  a  positive  lower  bound, 
further  efficiencies  can  be  gained,  in  com.bination  with  the 
procedures  of  Section  V.C  above.   Since  multiplication  is 

computationally  faster  than  division,  the  value  1/(A_/A  )  =  A  /A_ 

* 

can  be  precomputed  and  stored.   Thus  if  U.  <_  A_/A  , 

U.  ,  =  U.-A  /X    can  be  computed  as  the  next  thinning  variate. 
Note  that  no  intensity  function  calculation  is  required. 
However,  if  U.  >  VA  ,  the  thinning  method  proceeds  with 

the  next  step,  evaluating  the  intensity  function  at  the 

* 
point  T.  and  determining  whether  or  not  to  thin  the  point. 
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Now,  further  information  is  known  about  U . .   Specifically, 

either  U^  >  A(T.)/A  ,  in  which  case  T.  is  thinned,  or 

*  *    * 

X/X      <  U^  <_  A(T.)/A  .   In  either  case  U.  ,  can  be  computed 

by  "scaling  up"  U. . 

Thus,  the  algorithm  for  recycling  uniform,  random 

variates  for  thinning  is  as  follows: 

it  * 

1.  If  U^  £  X/X    ,  let  U^_^^  =  U^-A  /X    and  exit. 

2.  If  A/A*  <  Ui  £  A(T*)/A*,  let  U  ^_^^   =    (  A*  •  U^-A)  /  ( A  (T* ) -A ) 
and  exit. 

■k  k 

3.  Otherwise,  U.  >A(T.)/a  ,  let 
U^^^  =  (A*-U^-A(T*))/(A*-A(T*)) . 

By  precomputing  A  /A_,  this  recycling  procedure 
requires  only  one  multiplication  in  the  case  where 
U.  <_  A_/A  .   Otherwise  one  multiplication,  two  subtractions 
and  one  division  are  required.   In  either  case  the  recycling 
procedure  is  generally  faster  than  generating  uniform 
random  variates  from  a  random  number  generator,  even  when 
a  logical  IF  statement  is  added  to  check  for  extreme 
values  ("small  bits"). 

2 .   Recycling  Of  Exponential  Variates 

This  section  applies  only  where  the  intensity  func- 
tion is  exponential  polynomial.   Here  the  possibilities 
are  less  promising.   In  the  general  case  where  no  lower 


bound,  A_  is  used,  the  following  algorithm  would  apply: 

If  E.  >  In  A   -  In  A(T.),  let  E.^,  =  E.  -  In  A   +  In  A (T. 
1  —               1         1+1     1  1 

Otherwise 


Let  E^_^^  =  In  (A*  -  A(T*))/(A  -expl-E^)  -  A  (T^) 
where  E .  is  a  unit  exponential  random  variate. 
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*  * 

In  the  first  case,  E .  >^  In  A   -  In  A  (T.)/  a  time 

* 
savings  would  generally  be  realized  since  In  A   could  be 

computed  once  and  stored  and  In  A(T.)  is  simply  the  value 

*       *2 
of  the  polynomial,  i.e.  a^  +  a,T.  +  a^T.  ,  which  must  be 


computed  in  any  event.   In  the  second  case,  however,  the 
cure  is  truly  worse  than  the  illness.   It  is  faster  to 
simply  generate  another  exponential  variate,  assuming  they 
are  called  in  arrays. 

For  there  to  be  a  time  savings,  however,  it  must 
be  possible  to  make  a  reasonable  prediction  of  the  number 
of  exponentials  which  must  be  generated.   Otherwise  an 
excessive  number  of  calls  to  the  random  number  generation 
subroutine  may  destroy  the  gains  made  through  recycling. 

In  the  case  where  a  lower  bound,  ^,  for  the  inten- 
I   sity  function  exists  and  is  positive,  it  is  possible  to 
determine  the  expected  number  of  exponentials  which  must 
be  generated.   Variates  are  reused  if  they  are  greater  than 
In (A  /X) .      That  is,  if  E^  >  In (A  /X) ,    then  E^^^  =  E^  -  In ( X  /X) 

Otherwise,  a  new  (i.e.  non-recycled)  variate  is  used.   Thus 

* 
the  probability  of  not  recycling,  p,  is  A_/A   and  the  number 

of  variates  required  is  binomial  with  mean  n  p,  where  n 

is  the  number  of  points  to  be  thinned. 

Empirical  results  for  the  five  test  cases  considered 

are  shown  in  Appendix  B.   Using  the  calling  rule  of  expected 

number  plus  four  standard  deviations  for  generating  thinning 

exponentials  yielded  inconclusive  results  as  compared  with 


48 


the  procedure  in  which  exponential  thinning  variates  are 
generated  in  arrays  with  no  recycling.   As  expected,  for 
larger  N  ,  (Cases  II,  III  and  V),  recycling  provided  a  slight 
time  advantage  (seven  percentage  maximum)  while  for  small 
N  ,  (Cases  IV  and  VI)  recycling  was  slower.   In  case  VI, 
recycling  caused  run  time  to  be  approximately  five  percen- 
tage greater  than  that  without  recycling.   Using  a  calling 
rule  of  expected  number  plus  two  standard  deviations  reduced 
the  disadvantage  slightly  to  four  percent.   The  reason  that 
recycling  can  cause  longer  run  times  than  not  recycling  is 
that  an  additional  logical  IF  statement  is  required  for  the 
recycling  program.   Again,  when  exponentials  are  used  for 
thinning  and  the  mean  number  of  points  to  be  thinned  is 
on  the  order  of  two  or  three  hundred,  it  is  probably  not 
worth  the  effort  to  recycle.   If  several  thousand  "thinnings" 
are  required,  the  savings  may  indeed  be  worthwhile. 

Results  were  somewhat  surprising  in  the  general 
class  of  intensity  functions  for  which  uniform  variates  are 
used  for  thinning.   It  was  expected  that  a  significant  savings 
would  be  realized  since  the  uniform  variates  can  be  recycled 
in  all  situations.   In  fact,  test  runs  were  run  in  which 
only  one  uniform  variate  was  generated  for  the  entire  run 
with  recycling  used  throughout.   The  only  program  statement 
which  added  time  was  a  logical  IF  statement  to  preclude 
dividing  by  zero  and  to  diminish  the  probability  of  small 
bit  usage.   As  expected,  some  bias  was  experienced  in  the 
mean  and  an  unusually  high  variance  was  noted,  indicating  a 
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low  degree  of  "fidelity"  to  the  true  non-homogeneous  Poisson 
process  being  simulated.   Of  particular  interest  however, 
were  the  results  on  execution  time.   Since  this  setup 
!|   essentially  gives  the  lower  bound  on  execution  time  for 
recycling  of  uniform  variates.  it  was  expected  that  signi- 
ficant time  savings  would  be  realized  in  comparison  to  no 
recycling.   In  practice,  however,  the  savings  were  minimal, 
with  a  maximum  of  only  three  percent  savings.   This  is 
attributed  to  the  efficiency  of  the  LLRANDOM  package  in 
generating  uniform  variates  with  the  logical  IF  statement 
being  only  a  secondary  cause  to  longer  execution  time. 

Appendix  B  shows  results  for  both  exponential  and 
uniform  thinning  variates. 

The  key  point  to  keep  in  mind  here  is  that  the  above 
results  reflect  the  case  where  thinning  variates  are  called 
in  arrays,  i.e.  many  at  a  time.   Thus  the  comparison  is 
between  a  very  "fast"  variate  and  recycling.   As  discussed 
above,  calling  by  arrays  results  in  considerable  time  savings 
compared  to  one-at-a-time  generation  (up  to  fifty  times 
faster).   Thus,  in  the  case  where  a  slower  random  number 
generator  is  utilized  or  where  variates  are  called  one  at 
a  time,  use  of  the  lower  bound  may  indeed  result  in  con- 
siderable time  savings. 

F.   FINAL  PROGRAM 

A  general  program  was  developed  which  incorporated  the 
efficiency  considerations  discussed  above.   The  program  is 
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general  in  that  it  can  be  used  with  any  of  the  general 
class  of  intensity  functions,  whether  exponential  polynomial 
or  not.   The  program  is  essentially  four  programs,  each 
used  in  a  specific  case.   The  program  classifies  simulations 
into  the  four  classes  by  asking  two  questions: 

1.  Is  the  intensity  function  exponential  polynomial? 

2.  Does  the  intensity  function  have  a  positive  lower  bound? 
The  first  of  these  determines  whether  uniform  variates 

(general  case)  or  exponential  variates  (exponential  poly- 
nomial case)  are  used  for  thinning.   The  second  consideration 
merely  deletes  an  unnecessary  logical  IF  statement  in  the 
case  where  no  lower  bound  is  used. 

The  computer  program,  NHPP ,  is  listed  after  the  appen- 
dices, and  requires  a  user  supplied  subprogram  FUNCTION  FCN(T) 
to  compute  the  intensity  function  values  for  each  value  of 
t.   If  the  intensity  function  is  exponential  polynomial, 

only  the  exponent  portion  should  be  calculated,  i.e.  the 

2 
statement  FCN  =  a^  +  a-,t  +  a^t   (for  degree-two  polynomial) 

should  appear  in  the  subprogram.   Otherwise  the  entire 

intensity  function  should  be  evaluated. 

Empirical  results  showed  that  the  final  program, 
utilizing  the  efficiency  considerations  mentioned  in  this 
chapter,  resulted  in  a  program  which  ran  in  two- thirds  the 
time  of  the  basic  thinning  program. 

Please  see  Appendix  B  for  case-by-case  results. 
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VI.   RESULTS,  CONCLUSIONS  AND  RECOMMENDATIONS 

A.   GENERAL 

This  section  presents  the  results  of  comparison  of 
the  Poisson  decomposition  and  gap-statistics  algorithm 
with  three  variations  of  the  thinning  algorithm.   These 
variations  are  the  basic  thinning  algorithm,  the  modified 
thinning  algorithm  (final  program)  and  the  special  case 
one-at-a-time  algorithm.   Section  B  presents  the  performance 
of  each  of  the  algorithms  when  measured  by  the  two  measures 
of  effectiveness,  computational  speed  and  computer  memory 
requirements.   Section  C  examines  the  results  with  a  view 
toward  identifying  the  strong  and  weak  points  of  each 
algorithm.   Section  D  recommends  further  avenues  of  study. 

Again,  in  comparing  the  two  classes  of  algorithms,  one 
basic  distinction  must  be  kept  in  mind.   That  is  that  the 
Poisson  decomposition  and  gap-statistics  algorithm  as  imple- 
mented by  PATROW  [Ref.  3]  is  limited  to  a  special  class  of 
intensity  functions,  i.e.  exponential  polynomial  of  degree 
two  (or  less) .   Although  the  algorithm  could  be  adapted  to 
higher  order  polynomials  (by  further  bisection  of  inter- 
vals) ,  the  already  complex  programming  considerations  would 
grow  significantly.   In  contrast,  the  thinning  algorithm 
is  a  completely  general  method  which  is  analytically  valid 
for  any  functional  form  of  permissible  intensity  functions 
(positive  and  right  continuous) . 
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The  results  presented  here  are  necessarily  limited  to 
that  class  of  intensity  functions  for  which  both  algorithms 
can  be  compared,  i.e.  the  degree-two  exponential  polynomial 
class.   The  purpose  was  to  determine  the  relative  performance 
of  the  thinning  algorithm  on  this  piece  of  common  turf  with 
the  heretofore  champion,  Poisson  decomposition  and  gap- 
statistics  . 

The  basic  result  is  that  the  thinning  algorithm  is 
indeed  quite  competitive  with  the  Poisson  decomposition  and 
gap-statistics  algorithm  in  the  area  of  mutual  validity. 
This,  combined  with  its  ease  of  programming  and  ability  to 
generate  variates  from  any  intensity  function,  make  the 
thinning  algorithm  a  highly  attractive  tool  for  generating 
non- homogeneous  Poisson processes  of  any  type. 

One  shortcoming  of  the  thinning  program  was  revealed  by 
the  first  test  case  considered  (see  Section  IV. C).  This  is 
a  fast  rising  exponential  polynomial  which  rises  from  a 

value  of  five  to  almost  3300  over  the  interval  (0,100].   For 

* 

X      =  3294.47,  the  expected  number  of  points  in  the  bounding 

process  is  329,447  while  the  non-homogeneous  Poisson  process 
being  simulated  has  an  expected  number  of  only  1464  points. 
Thus  all  but  one  point  in  about  200  are  thinned  out.   The 
thinning  algorithm  could  be  more  efficiently  adapted  to 
this  case  by  partitioning  the  interval,  alleviating  the 
necessity  to  store  over  300,000  bounding  process  points. 
However,  the  efficiency  involved  would  still  be  low,  and  the 
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best  solution  appears  to  be  to  utilize  the  Poisson 
decomposition  and  gap-statistics  approach.   The  key  point 
Ij   here  is  that  the  problem  is  easily  recognized  beforehand, 
as  discussed  in  Section  V.A,  and  avoidable. 

Table  II  presents  a  general  comparison  of  the  two 
algorithms . 

B.   MEASURES  OF  EFFECTIVENESS  RESULTS 

Chapter  four  details  the  comparison  procedure  utilized 
to  develop  the  following  results. 

1.  The  Basic  Thinning  Algorithm 
Salient  features  for  this  case  include  the  use  of 

uniform  variates  for  thinning  and  one-at-a-time  generation 
of  exponential  and  uniform  variates. 

Table  III  presents  the  results  for  each  of  the  five 
test  cases  run.   Algorithm  A  is  computer  program  DEGTWO, 
the  Poisson  decomposition  and  gap-statistics  program 
developed  by  PATROW  and  listed  at  the  end  of  this  paper. 
Algorithm  B  is  computer  program  NHPTHN,  the  basic  thinning 
algorithm,  also  listed. 

The  thinning  algorithm  was  fastest  in  two  out  of 
the  five  test  cases  run  and  required  eighty  percent  of 
the  core  space  required  for  the  gap-statistics  algorithm. 

Table  VII  lists  core  storage  requirements  for  each 
algorithm. 

2 .  The  Modified  Thinning  Algorithm  (Final  Program) 
This  section  compares  the  best  case  performance 

for  both  algorithms. 
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The  modified  thinning  algorithm  includes : 

1.  Use  of  exponentials  for  thinning  of  exponential 
polynomial  intensity  functions 

2.  Use  of  lower  bounds 

3.  Partial  recycling  of  exponential  thinning  variates 

4.  Use  of  exponential  variates  left  over  from  generation 
of  the  bounding  process. 

Each  of  these  refinements  are  discussed  in  Section  V. 

The  Poisson  decomposition  and  gap-statistics 
algorithm  used  was  again  the  implementation  by  PATROW 
[Ref.  3],  program  DEGTWO.   In  addition  to  the  normal  running 

of  the  program,  a  second  set  of  comparisons  was  made  uti- 

* 

lizing  separately  calculated  values  for  c  /  the  bound  for 

the  conditioning-acceptance-rejection  routine.  This  is 
discussed  by  PATROW  [Ref.  3].  These  runs  are  indicated 
by  asterisk  (*)  in  Table  IV. 

In  the  first  case,  the  thinning  algorithm  was  faster 
in  four  out  of  the  five  test  cases,  with  the  best  relative 

performance  occurring  in  Case  VI. 

* 

^'Then  the  improved  values  of  c  were  incorporated 

for  Cases  IV,  V,  and  VI,  the  Poisson  decomposition  and  gap- 
statistics  algorithm  improved  substantially,  winning  in 
three  out  of  the  five  test  cases. 

The  relative  advantage  in  computational  speed  was 
less  than  a  factor  of  two  in  all  cases  with  a  maximum  of 
1.83  to  1  in  favor  of  the  thinning  algorithm  in  case  VI. 
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Core  storage  for  the  thinning  algorithm  was  determined 
by  using  only  the  part  of  the  algorithm  which  used  exponen- 
tial thinning  variates.   This  precluded  the  requirement  for 
storing  5000  uniform  variates  which  are  not  used.   The 
Poisson  decomposition  and  gap-statistics  program  (DEGTWO) 
required  about  88,000  bytes  of  core  storage  as  compared  to 
about  94,000  for  the  thinning  program  (NHPP  modified  to 
exclude  unused  uniform  variate  array) .   Detailed  results 
are  listed  in  Tables  IV  and  VI. 

3 .   The  One-At-A-Time  Thinning  Algorithm 

As  discussed  in  Section  II. C,  the  one-at-a-time 
algorithm  was  developed  only  to  test  the  relative  efficiency 
of  the  algorithm  used  to  generate  the  next  event  in  a  non- 
homogeneous  Poisson  process.   This  latter  requirement  may 
arise  in  event-step  simulation  where  only  the  next  event 
in  a  non-homogeneous  Poisson  process  is  desired  rather  than 
an  array  containing  all  events  in  a  specified  interval. 

Computationally,  the  one-at-a-time  algorithm  is 
quite  similar  to  the  basic  thinning  algorithm.   The  only 
essential  difference  is  that  the  basic  thinning  algorithm 
generates  and  stores  all  the  points  in  the  bounding  process 
(intensity  function  A  )  before  thinning,  whereas  the  one- 
at-a-time  algorithm  generates  a  point  in  the  bounding  process 
and  thins  that  point  before  continuing.   The  latter  method 
removes  the  requirem.ent  for  an  additional  array  to  store 
the  bounding  process  points .   This  in  turn  saves  about 
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20,000  bytes  of  core  storage  requirement  when  the  programs 
are  dimensioned  to  accept  5000  points. 

As  implemented  here,  the  one-at-a-time  algorithm 
sim.ply  generates  the  next  point  in  the  non-homogeneous  Poisson 
process  and  stores  it,  stopping  when  the  last  point  generated 
lies  outside  the  interval.   All  variates  in  this  program  are 
generated  one  at  a  time.   The  results  shown  in  Table  V  are 
thus  a  good  indicator  of  the  relative  efficiency  of  using 
this  method.   As  can  be  seen,  the  one-at-a-time  algorithm 
(program  NHPOAT)  is  faster  than  the  Poisson  decomposition 
and  gap-statistics  algorithm  (program  DEGTWO)  in  three  of 
the  five  test  cases  run.   This  is  true  despite  the  fact  that 
DEGTWO  generates  all  variates  in  arrays,  taking  advantage 
of  the  time  economies  of  scale  mentioned  in  Section  V.   The 
one-at-a-time  algorithm  also  requires  forty  percent  less 
core  space. 

From  the  tables  one  can  also  see  that  since  both  the 
best  case  thinning  algorithm  (Table  IV)  and  the  one-at-a- 
time  thinning  algorithm  (Table  V)  are  compared  to  the 
Poisson  decomposition  and  gap-statistics  algorithm,  it  is 
possible  to  obtain  a  reasonable  comparison  of  the  best  case 
thinning  algorithm  and  the  one-at-a-time  thinning  algorithm. 
For  example,  for  the  sample  intensity  function  used  in 
Case  II  (see  Figure  3),  the  ratio  (. 8127/ . 5541)  =  1.47 
indicates  that  execution  time  for  the  one-at-a-time  thinning 
algorithm  is  almost  fifty  percent  greater  than  that  of  the 
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best  case  thinning  algorithm.   The  ratios  for  Cases  III, 
IV,  V  and  VI  are  1.62,  1.27,  1.37,  and  1.21  respectively. 

The  tables  also  demonstrate  the  time-of-day  effect 
discussed  in  Section  IV.   That  is,  the  execution  time  for 
a  given  program  is  not  the  same  each  time  it  is  run.   For 
example,  program  DEGTWO  took  21.8808  seconds  for  the  run 
recorded  in  Table  V  compared  to  18.8917  seconds  for  the  run 
recorded  in  Table  IV.   For  this  reason,  ratios  of  execution 
times  were  chosen  as  the  measure  of  effectiveness  rather 
than  absolute  times. 

A  FORTRAN  program  listing,  NHPNXT,  is  provided  at 
the  end  of  this  thesis.   This  program  generates  the  next 
point  in  a  non-homogeneous  Poisson  process  with  a  user 
supplied  intensity  subprogram,  FUNCTION  FCN.   This  program 
can  be  used  in  conjunction  with  event-step  simulation 
programs,  including  SIMSCRIPT,  where  it  is  desired  to  mini- 
mize core  space  (at  the  expense  of  speed)  or  where  only  one 
event  is  desired  at  a  time.   Core  requirements  are  shown 
in  Table  VI. 

C.   CONCLUSIONS 

Both  the  thinning  algorithm  and  the  Poisson  decomposition 
and  gap-statistics  algorithm  include  two  general  types  of 
operations:   a  "generating"  process  and  a  "second  stage". 
For  the  Poisson  decomposition  and  gap-statistics  algorithm, 
the  generating  process  is  the  non-homogeneous  Poisson  process 
with  degree-one  exponential  polynomial  intensity  function. 


62 


<u 
-p 
>1 

m 

•H 


>1 

o 

e 

0) 

s 

+J 
:3 
a 
g 
o 
o 


03 

2 


en 
O 


•H 
O 

cn 

rH 


CN 

^ 

o 

•^ 

LD 

m 

iH 

o 

r- 

rH 

H 

H 

•. 

^ 

«k 

r- 

o 

'* 

O 

CO 

r~ 

rH 
rH 

ro 

^-> 

cn 

^-^ 

04 

CD 

z 

P^ 

■P 

ffi 

ffi 

(0 

1 

&H 

z 

^ 

a 

0^ 

— ' 

CJ 

(U 

X 

G 

U 

2 

e 

0) 

■ — ' 

x: 

o 

T3 

-p 

C 

e 

■H 

_r^ 

fC 

^ 

S-l 

u 

4J 

o 

H 

C 

•H 

D^ 

-G 

0 

S-l 

H 

S 

•H  .— ~ 

0 

< 

4J  O 

en 

e  Eh 

•H   S 

H 

m 

x:  X 

en  Eh 

< 

G 

-P  2 

0  U 

■H 

•H   CU 

a,a 

CJi 

C 

U   E 

B  Q 

c 

G 

0  2 

0  — 

•H 

■H 

Cn  — 

o 

G 

-C 

H 

QJ    to 

G 

Eh 

-=3:  ^ 

Q    U 

H 

G 

•H 

X: 

^ 

Cn-H 

C  -P 

Eh 

0) 

G    0 

0  tn 

H 

•H    Oi 

W  -H 

0 

m 

G 

ui  -u 

H 

H 

G   -P 

•H    03 

W 

T3 

•H    X 

0  -P 

td 

0 

x:  Q) 

a.  en 

ffl 

S 

Eh  Z 

H 
> 

EH 


CO 
Eh 

W 

s 

H 

D 

a 
w 

>H 

Pi 

o 

w 

W 
Eh 
D 

Oi 

s 
o 
u 


CM 


m 


63 


For  the  thinning  algorithm  (as  implemented  herein)  the 
generating  process  is  homogeneous  Poisson. 

The  second  stage  for  the  Poisson  decomposition  and  gap- 
statistics  algorithm  is  the  actual  decomposition  and  genera- 
tion of  variates  by  the  conditioning-acceptance-rejection 
method.   For  the  thinning  algorithm,  the  second  stage  con- 
sists of  the  thinning  of  the  points  in  the  bounding  process. 
Thus  one  algorithm  generates  events  and  adds  more  events 
from  a  second  process  while  the  other  generates  events  and 
subtracts  some  out. 

The  strongest  point  in  the  Poisson  decomposition  and 
gap-statistics  algorithm  is  the  highly  efficient  generation 
of  the  events  in  the  degree-one  exponential  polynomial  inten- 
sity function  process.   This  is  done  with  the  gap-statistics 
algorithm  which  is  two  to  five  times  faster  than  the  thinning 
algorithm  for  this  type  of  process  (see  Appendix  A) .   At 
the  same  time,  the  conditioning-acceptance-rejection  routine 
is  relatively  quite  inefficient. 

There  are  many  considerations  in  predicting  the  relative 
success  of  the  two  algorithms,  i.e.  which  will  be  faster  in 
a  given  situation.   For  example,  the  Poisson  decomposition 
and  gap-statistics  algorithm  is  affected  by  factors  such 
as  whether  or  not  partitioning  is  required,  the  percentage 
of  the  total  number  of  variates  which  come  from  the  degree- 
one  exponential  polynomial  process,  and  whether  time  rever- 
sal is  required.   For  the  thinning  algorithm,  the  fraction 
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* 

of  the  lower  bound  divided  by  the  upper  bound  \/\    ,    would 

seem  to  be  a  good  indicator  of  success. 

For  the  test  cases  considered,  however,  the  only  con- 
sistent indicator  was  the  expected  number  of  events  in  the  non- 
homogeneous  Poisson  process  being  simulated.  The  smaller  the  number 
of  events  in  the  non-homogeneous  Poisson  process  being  simu- 
lated, the  better  the  relative  performance  of  the  thinning 
algorithm  over  the  Poisson  decomposition  and  gap-statistics 
algorithm..   Thus  it  appears  that  each  algorithm  has  a  fixed 
and  variable  part  in  terms  of  time.   The  thinning  algorithm 
has  a  shorter  "setup"  cost  in  terms  of  time  but  the  variable 
cost  or  "cost  per  additional  variate"  seems  to  be  smaller 
for  the  Poisson  decomposition  and  gap-statistics  algorithm. 

The  exact  cause  of  this  phenomenon  is  not  known  although 
it  appears  to  be  centered  in  the  conditioning-acceptance- 
rejection  routine. 

In  the  larger  spectrum  of  non-homogeneous  Poisson  process 
generation,  it  seems  clear  that  the  thinning  algorithm  is 
the  best  all-around  method  available. 

D.   RECOMMENDATIONS 

Two  specific  areas  for  further  study  are  recommended. 

First,  the  thinning  algorithm  as  implemented  here  uses 
only  homogeneous  Poisson  processes  for  bounding.   It  might 
be  worthwhile  to  investigate  the  possibility  of  using  other 
processes,  such  as  non- homogeneous  Poisson  processes  with 
degree-one  exponential  polynomial  intensity  functions,  as 
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bounding  processes.   This  would  allow  the  efficient  gap- 
statistics  algorithm  to  be  utilized  although  function 
evaluation  would  in  general  become  more  time  consuming. 

The  second  area  is  in  finding  the  optimum  method  for 
generating  the  degree-two  exponential  polynomial  class  of 
intensity  function.   These  will  undoubtedly  remain  of 
interest  due  to  their  statistical  properties.   Here,  it 
seems  clear  that  the  best  features  of  the  two  algorithms 
can  be  combined.   Specifically,  the  Poisson  decomposition 
and  gap-statistics  algorithm  can  be  modified  to  use  thinning 
rather  than  conditioning-acceptance-rejection  for  generating 
the  points  in  the  difference  function  process,  i.e.  the 
process  with  intensity  function  A   =  A (t)  -  A  (t) .   Also 
the  criterion  for  the  decomposition  might  preferably  be 
that  the  intensity  function  of  the  remainder  be  monotonically 
increasing.   This  would  make  it  easy  to  find  the  upper 
bound  for  the  function,  and  the  most  efficient  version  of 
the  thinning  algorithm,  that  where  the  number  of  computations 
of  the  intensity  function  is  minimized,  could  be  used. 
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APPENDIX  A 

GENERATION  OF  DEGREE-ONE  EXPONENTIAL 
POLYNOMIAL  INTENSITY  FUNCTIONS 

The  generating  process  for  the  Poisson  decomposition 
and  gap-statistics  algorithm  is  a non- homogeneous  Poisson 
process  with  degree-one  exponential  polynomial  intensity 
function.  This  is  generated  by  using  the  gap-statistics 
method,  which  is  subroutine  NHPP2  in  the  DEGTWO  program 
(see  listing  below) . 

To  determine  the  relative  speed  of  the  thinning  algorithm 
compared  to  the  gap-statistics  algorithm,  two  simple  degree- 
one  exponential  polynomial  intensity  functions  were  developed 
and  simulated. 

Table  VII  presents  the  results.   Case  A  is  a  monotone 
decreasing  intensity  function,  A (t)  =  exp(3.4  -  0.02-t) 
over  the  interval  (0,100].   Case  B  is  a  monotone  increasing 
intensity  function,  .\(t)  =  exp(.693  +  0.03-t)  over  the 
interval  (0, 50]  . 

Results  show  that  the  gap-statistics  algorithm  is  from 
two  and  a  half  (Case  B)  to  four  and  a  half  (Case  A)  times 
faster  than  the  thinning  algorithm. 
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APPENDIX  B 
RESULTS  OF  EFFICIENCY  MODIFICATIONS 

This  appendix  presents,  in  tabular  form,  the  results 
of  comparison  of  the  programming  modifications  listed  in 
Section  V. 

Table  VIII  shows  the  effects  of  utilizing  lower  bounds 
for  the  intensity  function.   Test  conditions  include: 

1.  Use  of  uniform  thinning  variates 

2.  Recycling  of  thinning  variates 

3.  Use  of  arrays  of  variates 

Table  IX  shows  the  gains  realized  by  employing  exponen- 
tial thinning  variates  in  contrast  to  uniform  thinning 
variates  for  exponential  polynomial  intensity  functions. 
Test  conditions  include: 

1.  Use  of  arrays  of  variates 

2.  Recycling  of  thinning  variates 

3.  Use  of  lower  bounds  for  intensity  function 
Table  X  shows  the  results  of  recycling  versus  no 

recycling  where  uniform  variates  are  used  for  thinning 
while  Table  XI  shows  the  same  comparison  when  exponential 
variates  are  used  for  thinning.   For  both  cases,  test 
conditions  include: 

1.  Use  of  arrays  of  random  variates 

2.  Use  of  lower  bounds  for  intensity  function 
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Table  XII  presents  the  results  of  incorporating  all 
of  the  programming  improvements  into  program  NHPP.   The 
final  thinning  program,  NHPP,  is  compared  to  the  basic 
thinning  program  without  modifications,  NHPTHN.   The 
essential  differences  are: 


NHPP  (final  program) 
Arrays  of  variates  generated 


Exponential  variates  used 
for  thinning 

Lower  bound  of  intensity 
functions  utilized 

Thinning  variates  recycled 


NHPTHN  (basic  program) 

Variates  generated  one 
at  a  time 

Uniform  variates  used 
for  thinning 

Lower  bound  =  0.0 


No  recycling  used 
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C      SUBROUTINE    NHPP 

C 

C  SUBROUTINE    NHPP 

C 

C  PURPQ  SE 

C  SIMULATES    A    NON-HOMOGENEOUS    POISSGN    PROCESS 

C  WITH    INTENSITY    FUNCTION    FCN(X}    OVER    A    GIVEN 

C  INTERVAL  USING  THE  THINNING  ALGORITHM. 

C 

C  USAGE 

C  CALL  NHPP(IStEL»ER,U8,XMI.NI,NTYPE,N,I£R) 

C 

C 

C  DESCRIPTION    OF    PARAMETERS 

C  IS    -    RANDOM    NUMBER    SEED.       ANY    INTEGER    WITH 

C  NINE    OR    LESS    DIGITS. 

C  EL    -    LEFT    END    POINT    OF    INTERVAL 

C  ER    -    RIGHT    END    POirJT    OF    INTERVAL 

C  UB    -    UPPER    BCUND    CF    THE    INTENSITY    FUNCTION 

C  FCN(X)     OVER    THE     INTERVAL     (£L,ER).        THE 

C  CLOSER    UB    IS    11   THE   LUB,    THE    NORE 

C  EFFICIENT    THE     PROGRAM 

C  XMIN    -    LOWER    BOUND    OF    THE    INTENSITY    FUNCTION 

C  OVER    THE    INTERVAL.       THE    CLOSER    XMIN     IS    TO 

C  THE    GLB,    THE    MORE    EFFICIENT    THE    PROGRAM. 

C  NTYPE   -    1     IF    THE    INTENSITY    FUNCTION     IS 

C  EXPONENTIAL    oolymomiaL»     1»E.    OF 

C  THE    FORM     EXP ( A    +    A     *X     ♦    A2*X**2    +...) 

C  0    OTHERWISE 

C  N    -    THE    TOTAL    NUMBER    OF    EVENTS    IN    THE 

C  NON-HOMOGENEOUS    POISSON    PROCESS 

C  lER    -    ERROR    FLAG.       lER    HAS    FOLLOWING    MEANINGS: 

C  1.  ..ER    I  S   LESS    THAN    cL 

C  2.  ..UB     IS    NON-POSITIVE 

C  3.  ..XMIN    IS    NEGATIVE 

C  4.. .MORE  THAN  5000  EVENTS  REQUIRED  FOR 

C  BOUNDING  PROCESS:  STORAGE  CAPACITY 

C  EXCEEDED. 

C  5.. .XMIN  IS  GREATER  THAN  UB 

C 

C  COMMENTS 

C  CALLING    PROGRAM    mjST    HAVE    A    COMMON    REGION, 

C  DCNNA,    OF   DINENSiaN(5000) 

C 

C  EXAMPLE:       DIMENSION    T(5000) 

C  COMMON/DQNNA/T 

C 

C  TIMES    TO    EVENTS    WILL    BE    STORED    IN 

C  CELLS    T(l)    THROUGH    T(N) 

C 

C  THE     INTENSITY    FUNCTION     IS    USER    SUPPLIED. 

C  IF    THE     INTENSITY    FUNCTION    IS    NOT    EXPONENTIAL 

C  POLYNOMIAL,     I.E.     IF   N^YPE    =    0,     THE    ENTIRE 

C  INTENSITY    FUNCTION    IS    EVALUATED. 

C 

C  EXAMPLE:     FUNCTION    FCN(X) 

C  A    =    1.0 

C  FCN    =    3=^X 

C  RETURN 

C  END 

C 

C  IF    THE    INTENSITY    FUNCTION    IS     EXPONENTIAL 

C  POLYNOMIALt     I.E.     NTYPE    =    1,     ONLY     THE 

C  PCLYNOMIAL     IS    EVALUATED. 

C 

C  EXAMPLE:    FLiNCTION    FCN(X) 

C  A    =    1  .0 

C  Al    =    0.5 

C  A2    =    0.05 

C  FCN    =    A    +    A1*X    i-    A2*X**2 

C  RETURN 
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c 
c 
c 
c 


c 
c 
c 


c 
c 
c 
c 


c 
c 
c 

c 
c 
c 
c 
c 
c 


c 
c 
c 
c 
c 
c 
c 


c 
c 
c 
c 
c 
c 
c 


END 

programmer:       LCDR    JOHN    SCOTT    REDD,     LSN 
SEPTEfVBER    1978 

SUBROUTINE    NHPP    (I S , EL, ER, UB» XMIN ,NTYPE ,N , I ER  ) 

OiWENSiGM    U{5C00),     TIMES(5000),    TT(5C00),    EE(5000) 

CGVVGN  /DCNNA/  TT 

COMMON  ANNE/  TIMES, EE 

EXTERNAL  FCN 

CALL  OVFLOW 

INITIALIZE    VARIABLES 

lER    =    0 
2K    =     .0001 
ZL    =    lOE-6 
PCTMIN    =     .05 

CHECK    FOR    PARAMETER    ERRORS 

IF  (EL.GE  .ER)    lER    =    1 

IF  (UB.LE.0.0)     lEP     =    2 

IF  (XMIN.LT.O  .0)     IER  =  3 

IF  (UB-XMIN.LT.ZL)     IER=5 

IF  (lEP.NE.O)    GO    TO     14 

GENERATE    POINTS    IN    HOMOGENEOUS     PCISSCN    PROCESS 
WITh    RATE    =    U3 

CALL    HPP     (IS,EL,ER,UB,NTY?E,NSTAR,NL£FT,NEXP, lEP) 
IF     {  IER.EQ.4)    GO    TQ     14 

IS    INTENSITY    FUNCTION    OTHER    THAN    LOG    CLADRATIC? 

IF     (NTYPt.EQ.G)    GO    TO    9 

LOG    gU4CRATiC    INTENSITY    FUNCTION 

DOES    IT    HAVE    A    MINIMUM    CR    IS    MINIMUM    LESS    THAN 
PCTVN    CF    MAX? 

PCT    =    XMIN/UB 

IF    (PCT. LT. PCTMIN)     GO    TO    6 

USE    MINIMUM 

CCMPUTt    EXPECTED    NUMBER    OF    EXPONENTIALS    NEEDED,     TAKING 
INTO    ACCOUNT    REUSEING    OF    THOSE    IN    MI N 

P    =    l.O-PCT 
Q    =    PCT 

XNST^R    =    FLOAT(NSTAR ) 

NCHK    =    IFIX(XNSTAR*P+4.0*SQRT(XNSTAR*P*0 n 
NCHK    =    MI N0( NCHK, NSTAR) 
CALL    REORD     (NCHK, N EXP ,NLEFT ) 
NCALL    =    NCHK-NLEFT 
IF    (NCALL  .LE.O)    GO    TO    1 
CALL    SEXPQN     (  IS  ,  EE,  NCALL  ) 
1     K    =    0 
KK     =    1 
EB    =    EE(KK) 
UBLN    =    ALQGCUB) 
BND    =    ALOGtUB/XMIN) 
KKK    =    0 

I    CCUN^S    THE    HPP    EVENTS 

K    CGUNTS    THE    NHPP    EVENTS 

KK    COUNTS    THE    CURRENT    LNI T    EXPONENTIALS     (FOR    THI^NING) 

KKK    COUNTS    THE    TOTAL    NUMBER    OF    UNIT    EXPONENTIALS 

(FCR    THINNING) 
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c 
c 
c 

c 
c 

c 


c 
c 
c 


c 
c 
c 


c 

c 
c 
c 


c 
c 
c 
c 
c 
c 


c 
c 

c 
c 
c 


c 
c 


00    5    1=1,  NSTAR 

IF     (EB.LT.BND)    GO    TO    2 

K    =    K+1 

TT(K)     =    TIMESd) 

EB    =    E3-BND 

GO    TO    5 

2  VAL    =    -FCN(TIMES(I  )  )+UBLN 
IF    (EB.LT  .VAL)    GO    TO    3 

K    =    K+1 

TT(  K)     =    TIMESd  ) 

3  KK    =    KK+1 

CHECK  TO  SEE  IF  MORE  UNIT  EXPO.M  NEEDED 

IF  (Kk.LE.NCHK)  GO  TQ  4 

GENERATE    "^OR  E    UNIT    EXPONENTIALS    FOR    THINNING 

KKK    =    KKK+KK 

FN    =    P*FLGAT(NSTAR-KKK) 

NEB    =    MAX0(1,IPIX( PN+4.0*SgRT (PN*C) ) ) 

CALL    SEXPGN     (IS,E£,NEB) 

KK    =    1 

NCHK    =    NEB 

4  £8    =    EE(KK) 

5  CONTINUE 
N    =    K 

G3    TO    14 

LOG    QUaORATIC    WITH    NO    MINIMUM 

6  CONTINUE 

CALL    REORO     (N STAR ,N EXP  ,NL£F T) 
NCALL    =    NSTAR-NLEFT 
IF     (NCALL. LE.O)     GO    TO    7 
CALL    SEXPON    ( I S, EE , NC ALL ) 

SET    VARIABLE  S 

7  K    =    0 

UBLN    =    ALOG(UB) 

I     COUNTS    HPP    EVENTS 
K    COUNTS    NHPP    EVENTS 

DO     8    1=1, NSTAR 

VAL    =    -FCN(TIM£S(I)  > +UeLN 

IF     (EE( I) .LT.VAL)     GC    tq    8 

K    =    K  +  1 

TT  (K)    =    T  IMES  (I  ) 

8  CONTINUE 
N    =    K 

GG    TC    14 

INTENSITY    FUNCTION    IS    NOT    LOG    QUADRATIC 

DOES    IT    HAVE    A    MINIMUM    CR     IS    MINIMUM    LESS    THAN 
PCTMN    GF    MAX? 

9  PCT    =    XMIN/UB 

IF     (PCT.lT.PCTMIN)    go    to    12 

USE    MINL'-UM 

INITIALIZE    VARIABLES 

NUMF    =    NSTAR 

CALL    SRANO    ( IS,U,NUNIF ) 

K    =    0 

FF    =    1.0/UB 

I     COUNTS    HPP     EVENTS 


c 
c 


c 
c 
c 


c 
c 
c 


c 
c 
c 

c 


c 
c 

c 
c 

c 


c 
c 

c 


K    COUNTS    NHPP    EVENTS 

DO    11    1=1  ,NST/5R 

IF     (U(I), GT.PCT)    GO    TO    10 

K    =    K+1 

TT(K)     =    TIMES (I ) 

GO    TO    11 

10  VAL    =    FCN(TJMES(I ) ) «FF 

!!=     (U(I  ).GT.  VAL)     GO    TO    11 

K    =    K+1 

TT(K)     =    -^IMfS  (I) 

GO    TO    11 

11  CONTINUE 
N    =    K 

GO    TO    14 

NO    MINIMUM 

12  NUMF    =    NSTAR 

CALL    SRAND    (IS,U,NUNIF) 

K    =    0 

FF    =    1.0/UB 

DO    13    1  =  1, NSTAR 

VAL    =    FCN(T  IMESd  )  )  *FF 

IF     (U(I  ).GT.  VAL)     GO    TO    13 

ACCEPT    POINT 

K    =    K  +•  1 

TT(K)    =    TIMES(I) 


REJECT    POINT 

13  CONTINUE 
N    =    K 

14  RETURN 
END 

SUBROUTINE    HPP 

SUBROUTINE    HPP    GENERATES    PCINTS     IN    A    hCMQGENEOUS 
PQISSON    PROCESS    WITH     INTENSITY    FUNCTION    =    US 

SUBROUTINE    HPP(IS,FLt5R,UB,NTYPE,NSTAR,NLEFT,NEXP,IER) 
DIMENSIG-J    TIM£S(5000),     EE(500C) 
CC>'l>'GN    /ANNE/    TIMES, tE 

INITIALIZE    VARIABLES 

5XMEAN    =     1.0/LB 

NEXP    =    IFiX(UB*(ER-EL)+4.0*SQRT(UB*(ER-EL) ) ) 

IF     (NEXP.GT,5000)     NEXP=5000 

CALL    SEXPON    (IS,EE,NEXP) 

SUM    =    EL 

ISTART    =    0 

ISTOP   =    NEXP 

00    1    JJ=1,NEXP 

E    =    EXMEAN*££(JJ  ) 

SUM    =    SUM+E 

TiME5(JJ)     =     SLM 

IF     (SUM.LT.ER  )    GO    TO    1 

NSTAR    =    JJ-1 

GO    TO    5 

1  CONTINUE 

EXCEPTIONAL     NUMBER    OF    POINTS    NEEDED 

2  ISTART    =    ISTART+JJ 

NEXP    =    IFIX(UB*(ER-SUM )+4.0*SQRT(UB*(ER-SUM) ) ) 
IF     (NEXP.EQWOJ    NEXP=1 
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iSTCP    =    ISTART^MEXP 

IF     (1ST0P.GT.5000)     GO   TO   4 

CALL    SEXPQN     (IStEE,NEXP) 

Ca    3    JJ=I ,NEXP 

E    =    eXMEAN*EE (JJ) 

SUM    =    SUM+E 

KK    =    JJ+I START 

TIMESiKK)     =    SUM 

IF     (SUM.LT.ER)    GO    TO    2 

NSTAR    =    KK-1 

G3    TQ    5 

3  CQNTIMUE 
C 

GO    TO    2 
C 

C   MORE  THAN  5000  EXPONENTIALS  NEEDED 
C 

4  i£R  =  4 
GJ  TO  6 

C 

C       CALCULATE    NUMBER    OF    EXPONENTIALS    NOT    USED 

C 

c 


5  NLEFT    =    NEXP-JJ 

6  RETURN 
ENC 

C. ..SUBROUTINE    REORD 

C         SUBROUTINE    REORD    REORDERS     EXPONENTIALS    LEFT    OVER    FROM 

C         SUBROUTINE    HPf>    FOR    USE    AS    THINNING    VAPIATES     IN    NHPP 

C 

SUBROUTINE    RECRD     ( NCHK, NEXP , NLEFT  ) 

DIMENSION    EE(=000),    TIVES(5000) 

CCMVON    /ANNE/    TIMES, EE 
C 

C      ARE    ENOUGH    EXPONENTIALS    LEFT    OVER    FRG^'    hPP? 
C 

IF     (NLEFT. GE. NCHK)     GO    TO    4 
C 

C       MCRE    EXPONENTIALS   NEEDED 
C 

C      REORDER    TOP    TO    BOTTOM    OR    VICE    VERSA? 
C 

IF     (NEXP. LT. NCHK)     GO   TG    2 
C 

C       RECROEP    BOTTOM    TQ    TOP 
C 

Nl    =    NCHK -NLEFT 

N2    =    NEXP-NLEFT 
C 

CO    1    1=1, NLEFT 

J    =    Nl+I 

K    =    N2+I 

E£(J)    =    EE(K) 

1  CONTINUE 
GO    TO    6 

C 

C      REORDER    TOP    TO    BGTTOM 

C 

2  NCHKO  =  NCHK+1 
NEXPO  =  NEXP-^1 
CO  2  1=1,  NLEFT 
J    =    NCHKO-I 

K    =    NEXPO-I 
EE( J )    =    EE(K) 

3  CONTINUE 
GO    TO    6 

C 

C       ENOUGH    LEFT    OVER    FROM    hPP    FOR    ALL   THINNING 


C 


Nl    =    NEXP-NLEFT 
DO    5    1=1, NCHK 
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J    =    Ni+1 

EE(  I)    =    EE(J) 

5  CONTINUE 

6  RETURN 
END 
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c 
c 

C  SUBROUTINE    NHPThM 

C 

C  PURPOSE 

C  SIMULATES     A    NQN  HOMOGENEOUS    POISSON    PROCESS 

C  WITH    INTENSITY    FUNCTION    FCN(X)    USING    . 

C  THE    THINNING    ALGORITHM 

C 

C  CALL    NHPTHN {  IS,ELt ER,U3,N,IER ) 

C 

C  DESCRIPTION    OF    PARAMETERS 

C  iS    -    RANDOM    NUMBER    SEED,      ANY     INTEGER    WITH    NINE 

C  OR    LESS    DIGITS. 

C  EL    -    LEFT    END    POINT    OF    INTERVAL 

C  ER    -    RIGHT    END    POINT    HF    THE    INTERVAL 

C  UB    -    UPPER    BOUND    OF    THE    INTENSITY    FLNCTION , FCN ( X ) , 

C  OVER    THE    INTERVAL     (EL,ER).       THE   CLOSER    UB    IS 

C  TO    THE    LEAST    UPPER    BOUND, LUB,    THE    MORE 

C  EFFICIENT    THE    PROGRAM.        UB    MUST    BE    STRICTLY 

C  POSITIVE 

C  N      -    THE    TOTAL    NUMBER    OF    EVENTS    IN    THE 

C  NON-HOMOGENEOUS    PCISSON    PROCESS. 

C  iER    -    ERROR    FLAG.        lER    HAS    FOLLOWING    MEANINGS: 

C  1...ER    IS    LESS    THAN    EL 

C  2...UB    IS    NON-=>OSITIVE 

C 

C         COMMENTS 

C  CALLING    PROGRAM    MUST    HAVE    A    COMMON    REGION, 

C  SCOTT,    OF    DIMENSI0N(5000) 

C 

C        EXAMPLE;   DIMENSION  T  (5  JOO  ) 

C  CCMMON/SCGTT/T 

C 

C  TIMES    TO    EVENTS    WILL    BE    STORED    IN    CELLS 

C  T(  I)    THROUGH    T(  M)  . 

C 

C  PROGRAMMER:       LCDR    JOHN    SCOTT    REDD,    USN 

C  AUGUST    1=73 

C 

SUBROUTINE  NHFTHN  ( I S , EL, £R,U B, N, I EP ) 

DIMENSION)  TIMES(5000),  TTT(50Q0) 

COf'NCN    /SCOTT/    TTT 

EXTERNAL    FCN 

CALL    QVFLOW 
C  INITIALI Z£    VARIABLES 

IER    =    0 

IF     (EL.GE.ER)     IER     =    1 

IF     (U3.LE.0.0)    IER    =    2 

1=     (lER.NE.O)     RETURN 
C 

C         GENERATE    POINTS    IN    HOMOGENEOUS    PCISSCM    PROCESS    WITH    RATE 
C         =    UB  . 
C 

EXM5AN    =    1.0/LB 

I    =    1 

SUM    =    EL 

1  CONTINUE 

CALL    EXPON     (IS,£, 1) 

E    =    E^EXMEAN 

SUM    =    SUM+E 

IF     (SUM.GT.ER  )    GO    TO    2 

TI  MESdl     =    SUy 

I    =    I  +  l 

GO    TO    1 

2  CONTINUE 
NSTAR    =    I-l 

C 

C         COMMENCE    THINNING    THE    NSTAR    POINTS 


C 


Fl    =   l.G/UB 
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K    =    0 

IF     (NSTAR.EG.O)    GO    TQ    4 

DO    2    I=1,NSTAR 

CALL    SRAND    (IStUa  I 

RATIO   =    FCNdlMESd  )  )*FI 

IF     (U.GT. RATIO)    GO    TQ    3 

K    =    K'4't. 

TTT(K)    =    TIMES(I) 

CONTINUE 

N    =    K 
RETURN 
N    =    C 
RETURN 
END 


I 
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C         SUBRCUTINE    MHPO/iT 

C 

C 

C      SUEROUTINE   NhPOAT 

C 

C  PURPOSE 

C  SIVULATES    A    NON-HOVQGENEOUS    PQISSGN    PROCESS 

C  WITH    INTENSITY    FUNCTION    FCN(X)     OVER    THE 

C  INTERVAL     (ELtER)     USING   THE    QNE-A T-A-TI ME 

C  THINNING    ALGORITHM 

C 

C... USAGE 

C  CALL    NHPOATdSrEL,  ER,UB,N,1ER) 

C 

C. ..DESCRIPTION  OF  PARAMETERS 

C 

C  IS    -    RANDOM    NUMBER    SEED.       ANY    INTEGER    WITH    NINE 

C  CR    LESS    DIGITS  . 

C  EL    -    LEFT    END    POINT    QF    INTERVAL 

C  ER    -    RIGHT    END    POINT    OF    THE    INTERVAL 

C  UB    -    UPPER    BOUND    OF    THE    INTENSITY     FUNCT  ION,  FC  N(  X  ) , 

C  OVER    THE    INTERVAL     (EL,ER).       THE    CLOSER    UE    IS 

C  TO    THE    LEAST    UPPER    BOUND, LUB,     THE    MORE 

C  EFFICl.ENT    THE     PROGRAM.       UB    MUST    BE    STRICTLY 

C  POSITIVE 

C  N       -    THE    TOTAL    NUMBER    OF    EVENTS    IN    THE 

C  NON-HO.MGGENcOUS    PCISSON    PROCESS. 

C  lER    -    ERROR    FLAG.        lER    HAS    FOLLOWING    MEANINGS; 

C  l..,£R    IS    LESS    THAN    EL 

C  2...UB     IS    NCN-POSITIVE 

C... COMMENTS 

C      CALLING  PRCGRAM  MUST  HAVE  A  COMMON  REGION,  DONNA,  OF 

C        DIMENSION( 5C00) 

C 

C  EXAMPLE:       DIMENSION    TT  (5000) 

C  CQMMQN/DONNA/TT 

C 

C  TIMES    TO    EVENTS    WILL    BE    STORED    IN    CELLS 

C  T(l)    THROUGH    T(N) 

C 

SUBROUTINE  NHPGAT  (  I  S  ,  EL,  E'R  ,U  B  ,N  ,  I  ER  ) 

DIMENSION  TT(5000) 

CCf'NCN    /OCNNA/    TT 

EXTERNAL    FC N 

CALL    CVFLOW 
C... INITIALIZE    VARIABLES 

lER    =    C 

IF     (EL.GE.ER)     lER    =    1 

IF     f UB.LE.0.0)     lER    =    2 

IF    (lER.NE.O)    RETURN 
C 

i    =    1 

EXMEAN    =     1.0/UB 

SUN    =    EL 
C 

C         GENERATE    POINT     IN    BOUNDING    PROCESS 
C 

1    CONTINUE 

COLL    EXPON     (  IS,E,1  ) 

E    =    E^EXMEAN 

SUM    =    SUM+E 

IF     <SUM.GT.ER)    GO    TG    2 
C 
C         THIN    THE    POINT 


C 


CALL    RANDOM     (IS,U,1) 

RATIO    =       FCN(  SUM  )*EXMEAN 

IF     (U.GT.RATIC)    GO    TO    1 

TT( I)    =    SUM 

1=1+1 

GO    TO    1 

N    =    I-l 
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RETLRN 
END 


C  SLBROUTINE    NHPNXT 

C 

C  SUBROUTINE    NHPNXT    GENERATES    THE    TIME    OF    THE    NEXT 

C  EVENT    IN    A    NCN-HOMOGENEGUS    POISSON    PROCESS    WITH 

C  RATE    FUNCTION    FCN(X)     (USER    SUPPLIED). 

C 

C  LSAGE: 

C  CALL    NHPNXT{IS,U6tGL8,XLAST,ER,XNEXT,IER) 

C 

C  DESCRIPTION    OF    PARAMETERS: 

C  IS           -    RANDOM    NUMBER    SEED.       ANY    INTEGER    WITH 

C  NINE    OR    LESS    DIGITS. 

C  UB           -    UPPER    BOUND    OF    THE    INTENSITY    FUNCTION 

C  OVER    THE    INTERVAL    (XLAST,ER». 

C  GL  E         -    GREATEST    LOWER    BOUND    OF    TJ-E    INTENSITY 

C  ^UNCTION    OVER    THE    INTERVAL    (XLAST,ER). 

C  SET    =    0    IF     UNKNOWN. 

C  XLAST    -    THE    TIME    OF   THE    LAST    EVENT    IN    THE    PROCESS 

C  ER            -    RIGHT    END    POINT    OF    THE    INTERVAL 

C  XNEXT    -    THE    NEXT    POINT    IN    THE    PROCESS.       IF    THERE    ARE 

C  NC    MCRE    POINTS    IN    THE    INTERVAL    (XLAST,ER), 

C  XNEXT    IS    ASSIGNED    THE    VALUE    ER    +    1.0    AND 

C  lER    IS    SET    AT     5. 

C  lER         -    ERROR    FLAG.       I ER    HAS    THE    FCLLOWING    MEANINGS: 

C  1...L3    IS    NON-POSITIVE 

C  2...XLAST    IS    GREATER    THAN    ER 

C  3.  ..GLB    IS    NEGATIVE 

C  5...  XNEXT    IS    GREATER     THAN    ER 

C 

C  COMMENTS 

C  THE    INTENSITY    FUNCTION,    FCN,     IS    USER    SUPPLIED, 

C  EXAMPLE:       FUNCTION    FCN (X ) 

C  FCN    =    1.0    ^    EXP  (-X) 

C  RETURN 

C  .END 

C 

c 

C  PROGRAMMER:       LCDR    JOHN    SCOTT    REDD,     USN 

C  AUG   1978 

C 

c 

C 

SUBROUTINE    NHFNXT    { I S, U B ,GLB, XLAST , ER,XNEXT , I ER ) 

EXTERNAL     FCN 

CALL    CVFLCW 

DATA    EXPO/0. 0/,U/0.0/ 

RMIN    =    GLB/UB 

lER    =    0 

IF     (UB.LE.0.0)    lER     =    1 

IF    (XlAST  .GE.ER)    I ER  =  2 

1=     (GL3.LT.0.C)     IER=3 

IF     (lER.EC.l.OR.  lER  .EQ.2)    RETURN 

IF     (IER.EG.3)    GLB    =    0.0 
C 

C         GENERATE    E*(I,J);       CHECK    TO    SEE    IF    ADDITION    OF    E*(I,J) 
C  EXCEEDS    ER 

C 

EXMEAN    =    1.0/UB 

XNE^«    =    XLAST 
C 

C    GENERATE  CNE  EXPONENTIAL  AND  SCALE 
C 

i  CALL  EXPON  ilS,  EXP0,1) 

EXPG    =    EXPQ^EXMEAN 

XNEW    =    XNEW+EXPO 

IF     (XNEW.GT.ER  )    GO    TO    3 
C 

C         GENERATE   UNIFORNfO,  1)    THINNING    VARIATE 
C 

CALL    RANDOM    (IS,U,1 ) 
C 
C         TEST    LOWER    BOUND    FOR    THINNING 
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c 
c 
c 


c 
c 
c 


IF    (U.LE.RMIN)    GO    TO    2 

TEST    =0R    THINNING 

RATIO    =    FCN{XNEW)/Ue 
1=    (U.LE. RATIO)    GO    TO    2 

GO    TO    1 

ARRIVAL    HERE    INCICATES    SUCCESSFUL    THIN.NING 

2  XNEXT    =    XNEVi 

RETURN 

ARRIVAL    HERE    INDICATES    NO    MORE    POINTS     IN    INTERV/SL 

3  lER    =    5 

XNEXT    =   ER>1.  C 

RETURN 

END 
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C  SLBROUTINE    DEGTWO 

C 

C  SUBROUTINE  DEGTWO 

C 

C  PURFCSE 

C  SIMULATES    A    NON-HOMOGENEOUS    POISSGN    PROCESS    WITH 

C  GUAORATIC    EXPONENTIAL     INTENSITY    FUNCTION    OVER 

C  A    GIVEN    INTERVAL    USING   THE    POISSC^-DECOMPOS  IT  ION 

C  AND    GAP     STATISTIC    ALGORITHM. 

C 

C  USAGE 

C  CALL    DcGTWO(IS,A, Al,A2tEL»ER,II ,N,IER) 

C 

C  OESCRIPTICN    OF    PARANETEFS 

C  IS    -    RANDOM    NUMBER    SEED.       ANY    INTEGER     WITH    NINE 

C  OP    LESS    DIGITS. 

C  A    -    CONSTANT    IN    INTENSITY    FUNCTION 

C  n    -    1ST    DEGREE    COEFF    IN    INTENSITY   FUNCTION. 

C  A2   -    2ND    DEGREE    COEFF    IN    INTENSITY    FUNCTION. 

C  EL   -    LEFT    END    POINT    OF    INTERVAL. 

C  £R    -    RIGHT    END    POINT    OF    INTERVAL 

C  II    -    0    FOR    TIMES    OF    EVENTS. 

C  1    FOR    TIMES    BETWEEN    EVENTS. 

C  N      -    A    VECTOR    OF    LENGTH    5.       N(l)    THROUGH    N(4) 

C  CONTAIN    NUMBERS    OF    EVENTS    FRCN    VARIOUS 

C  COMPONENTS    OF    THE    DECOMPOSED    INTENSITY 

C  FUNCTION.       N(5)     CONTAINS    THE    TOTAL    NUMBER 

C  OF    EVENTS    IN    THE    NON-HOMOGENEOUS    POISSCN 

C  PROCESS. 

C 

C  COMMENTS 

C  CALLING  PRCGRAM  MUST  HAVE  A  COMMON  REGION,  HOLD, 

C  OF  DIMENSION  (50C0),ANO  AN  INTEGER  ARRAY  OF 

C  DIMENSION  (5). 

C 

C  EXAMPLE:       DIMENSION    T(5O0O),N(5) 

C  COMMON/HOLO/T 

C 

C  CALLING  PROGRAM  MUST  CONTAIN  THE  FOLLWOING 

C  ASSIGNMENT  STATEMENT: 

C 

C  M=N(5) 

C 

C  CALLING  PROGRAM  MUST  USE  THE  FCLLCWING  J CL  CARDS 

C 

C  //       EXEC       FORTCLGt  IMSL    =    DP 

C  //FORT.SYSIN       DO      * 

C 

C  TIMES    TO    EVENTS    OR    TIMES    BETWEEN    EVENTS    WILL    BE 

C  STORED    IN    CELLS    T(l)     THROUGH    T(M). 

C 

C 

SUBROUTINE    DEGTWO     ( I S  , A  ,A1 , A2 t EL ,E R , 1 1 , N* I ER ) 

CiMENSiON    T'IMES(5000  )♦     T(5000),    N(5),    d(5) 

COMMON    /MIKE/    TIMES/HOLO/T 

CALL    OVFLOW 
C 

C  INITIALIZE    VARIABLES 


C 


C 
C 


P(l  )  =  A 

P(2)  =  Al 

F(  3  )  =  A2 

F(4)  =  0. 

P(  5)  =  0, 


DO     i    1=1, 5 
1    N(I  )    =    0 


C 

C 

C  IF    RATE    FUNCTION     IS    LESS    THAN    DEGREE     TV^O 

C  USE    NHPP2    ROUTINE    ONLY 
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c 

c 


c 


c 
c 


c 


If    (A2.EQ.0.)    GO    TO    2 

GO    TO    4 

CALL    NHPP2     (IS, a,ER,A,Al,Il,Nl,IER) 

N(5)     =    Nl 

IF    (Ni.EQ.O)     RETURN 


DO     3    1=1, Nl 
TIMESd)    =    Td) 
3    CONTINUE 


RETURN 
C 

C         DETERMINE   COEFFICIENTS    FOR    MODIFIED 
C         DEGREE    ONE    RATE    FUNCTION 
C 

4  TEST    =    -Al/  (2.*A2) 
TINT    =    ER-EL 

IF     (Al.GE,0.,AND.A2.GT.O. )     GO    TO    5 
GO    TO    6 

5  8=    A-A2*TINT**2 
Bl    =    Al+2.=i«A2*TINJT 
GO    TO    10 

6  B    =    A 

IF    ({Al  .LE.O.  .AND.A2.LT.0. ) .OR. ( Al . GT. 0. . AND . A2 .LT. 0. 
ITINT) )    GO   TO    7 
GO    TO    8 

7  Bi    =    A1-«'A2*TINT 
GO    TG   10 

a    IF     (A1,GT«0. . AND. A2.LT.0. .AND. TEST. LT.TINTi    GC    TO    9 
El     =    Al 
GO    TO    10 
9    81    =    Al/2. 
C 

C         MUST    THE    INTERVAL    BE    PARTITIONED? 
C 

10  IF     (A1*A2. LT.C.«AND. TEST. LT. TINT)     GOTO    11 
ERNEW    =    E  R 

GO    TO    12 

11  ERNEW    =    TEST^-EL 
C 

C         GENERATE    DEGREE    ONE    NHPP    ON    INTERVAL 


12  BB  =  B 
831  =  81 

CALL  NHPP2  (iS,EL,ERNEW,BB,3Bl,0,Nl ,IER) 
N(l )  =  Nl 
IF  {N(l).EQi.O  J  GO  TO  14 


DO  13  1=1, Nl 
TIM£S( I)  =  T(I  ) 
13  CONTINUE 
C 

C 

C         COMPUTE    LENGTH    OF    INTERVAL    AND    DETERMINE   VALUE 

C         OF    CSTAR    FOR    USE    IN    REJECTION    ROUTINE 


14    Q    =    ERNEW-EL 
El    =    A 

E2    =    A2*Q**2 
E3    =    A1*Q 

E4    =    Ai**2/ (4.*A2) 
E5    =    Al**2/(i2.*A2  ) 

IF     {Al.GE  .0. .AND.A2.GT. 0. )     GO    TO    15 

IF  (Al.LT  .0  ..AND. A2  .GT.O.  .ANO.TEST  .C-E  .TINT  )  GC  TO  16 
IF  (Al.LT.O.. AND. A2.GT.0.. AND. TEST.LT. TINT)  GC  TC  17 
IF    (Al.LE .O..AND. A2.LT.0. )     GO    TO    18 

IF     (A1.GT.0..AND.A2.LT.0..AND.TEST.GE.TINT  )    GO    TO    19 
CSTAR    =    E XP(E1-E4)-EXP( El) 
GO    TO    20 
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15  C5TAR    =    EXP(;ei)-5XP(El-E2) 

GO    TO    20 

16  CSTAR    =    EXP(EH-E2  +  E3)-EXP(E1  +  E3) 
GO    TO    20 

17  CST4R    =     EXP(E1-E4)-EXP(E1-E5) 
GO    TO    20 

18  CSTAR    =    EXP(E1»-EXP(E1+E3+E2) 
GO    TC    20 

19  CSTAR    =    EXP(E1+E3+E2)-EXP(E1) 
C 

C         COMPUTE    INTEGRAL    OF    MCCIFIED    DEGREE    TWO    RATE    FUNCTION 

C         OVER     INTERVAL 

C 

20  CALL    HELP    _(A,AliA2,eL,ERNEWTPMTR) 

C 

C    IDENTIFY  AS  FIRST  SU6INTERVAL 

C 

NOTE  =  1 
C 

C         GENERATE    REAL-IZATION    ON    POISON     (PMTR)     VARIATE 
C 


C 

c 


c 


PM 


LL    HELP     (A,Ai  ,A27  bLtiRNEWTPMTR) 

TR  =  PMTR-(EXP(B)«(EXP(B1*£RNEW)-EXP(B1*EL)  )  )/Bi 


21  CALL  PVAR  (IStP^TR,M) 
IF  (NCTE.EQ.l  )  GO  TO  22 
GO  TO  25 

C 

C  REJECTION    ROUTINE   USED    ON    FIRST    SU8INTERVAL 

C 

22  N{2)  =  M 
P(4)  =  B 
P{ 51    =    Bi 

IF     (N(2  ).EQ<.0  )    GO    TO    24 

CALL  REJECT  ( IS ,E L ,CSTA R  ,  P, Q, N (2)  ) 


DO  23  1=1, M 
TIMES(N(1 )+!)  =  T( I ) 

23  CONTINUE 
C 

C 

C    HAS  THE  INTERVAL  BEEN  PARTITIONED? 

C 

24  IF  (ERNEW.EQ.ER)  GO  TO  34 
GO  TO  27 

C 

C         USE    REJECTION    ROUTINE    ON    SECOND    PART     CF    INTERVAL 

C 

25  N(4)  =  M 
P(  4)  =  B 
P(5  )    =    81 

C 

C         IF    NO    EVENTS    OCCURRED    BYPASS    REJECTION    ROUTINE 

C 

IF     (N<4)  .EQ-.O  )    GO    TO    35 

Q    =    ER-ELNEW 

CALL    REJECT     {  IS , ELNEW, CST ARi P , Q , N ( 4 U 
C 

C         COPY    TI^^ES    OF    EVENTS    INTO    'TIMES'    ARRAY 
C 


N4    =    N(1)+N(2)-^N(3) 


DO     26    1=1, M 

TIN£S(N4-»-I)    =     T(  I) 
26    CONTINUE 
C 
C 

C         GENERATION    OF    VARIATES    COMPLETE. 
C         GO    TO    ORDERING    ROUTINE 


GO    TO    35 
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C    INTERVAL  PARTITION  WAS  REQUIRED.   MUST  NOW 

C    CONSIDER  SECOND  SUBINTERVAL 

C 

C         DETERMINE    COEFFICIENTS    FOR    MODIFIED 

C         JEGREE    CNE    RATE    FUNCTION 

C 

27  IP     (Al.GT.O. .AND.A2.lt. 0. )     GO    TO    28 
e    =    A-A2*TINT=**2 

Bl    =    AH-2.=«=A2*TINT 
GO    TO    2  9 

28  8    =    A+(A1  /2.  )*TINT 
81    =    Al/2.+A2*TINT 

29  ELNEW    =    ERNEW 
C 

C         GENERATE    DEGREE    ONE    NHPP    ON    INTERVAL 
C 

BB^  8 

BBi  =  Bl 

CALL  NhPP2  (IIS, ELNEW, £R, BB,  3Bl,0,N3tIER) 

N(5)  =  N3 

IF  (N(3).£0.0)  GO  TO  51 

N3  =  N(l)+N(2  ) 
C 

C         TR/JNS1^£R    TIMES     BETWEEN    ARRAYS 
C 
C 


00    30    1=1, N3 

TIMES(N3^I)     =    T(I  ) 
50    CONTINUE 
C 
C 

31  Q    =    TINT 
C 

C         DETERMINE    VALUE    OF    CSTAR     FOR    USE    IN 

C         THE    REJECTION    ROUTINE 

C 

E2    =    A2*Q**2 

E3     =    Al*q 

IF     (Ai.GT.O.,AND. A2.LT.0. )     GO    TO    32 

CSTAR    =    EXD(E1-E4)-EXP(E1-E5-E3-E2) 

GO    TC   33 

32  CSTAR    =    £XP(  El-E4)-EXP{Ei<-E3-«-E2) 
C 

C         COMPUTE    INTEGRAL    CF    MCCIFIED    DEGREE    TWG    RATE 

C         FUNCTION    OVER    SECOND    INTERVAL 

C 

33  CALL    HELP     ( A, Al , A2 r ELNEW ,ER , PMTR ) 

PMTR    =    PMTR-( EXP(B )*(£XP(31^£R)-EXP{81*ELNEW)  )  )/ei 
C 

C       IDENTIFY    AS    SECOND    SUBINTERVAL 
C 

NOTE    =   2 

GO    TO    21 
C 

C    PARTITION  OF  INTERVAL  NOT  REQUIRED,   COMPUTE  TOTAL 
C    EVENTS  AND  SUPERPOSE  TWO  EVENT  STREAMS 
C 

34  N(5  )  =  N(1)-HN(2) 

IF  (N(2)  .EQ.O)  GO  TO  38 

L3GN  =  N(l).+1 

JBGN  =  1 

CALL  COL  ATE  (  L3GN  ,N  (  5  ),  1 ) 

GO  TO  38 
C 

C         PARTITION   WAS    REQUIRED.       DETERMINE 
C      AMOLNT    OF    SORTING    NEEDED 

35  N(5)    =    N(l )+N (2 )+N(3 )+N(4) 

IF     IN(2). EQ.O. AND. N(4) • EQ.O)    GO    TC    38 
IF    (N(4).EQ  .0)    GO    TO    36 
If     {N(2  ).  EQiOl    GO    TC    37 
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C         MUST    SUPtRPQS-E    FOUR    EVENT    STREAMS 
C 

L3GN    =    N(l)-«-l 

Lf  I^    =    N(l  )4-N(2) 

CALL    COLATE    ( LBGN , LFIN» 1 ) 

L3GN    =    LF  IN>N (3)+l 

J3GN    =    LFiN^l 

CALL    COLATE     (  LBGN  ,N  (  5)  ,  JBGN  ) 

GG    TO    38 
C 

C         MIST    SUPERPOSE    FIRST    HALF    OF    ARRAY    ONLY 
C 

36    N2    =    N(l)+N(2) 

L3GN    =    N(  1)  +  1 

CALL    COLATE    (LBGN»N2fl) 

GO    TO    58 
C 

C  MLST    SUPERPOSE    SECOND    hALF    OF    ARRAY    ONLY 

C 

3  7    KK    =    N  ( 1 )  ♦N  ( 2  )  +1 

LBGN    =    N(l^i-N(2)^■N(3)^■1 

LFIN    =    N(  5) 

CALL    COLATE     (  LBGN , N ( 5  J ,  KK ) 

GO    TO    5  8 
C 

C         ARE    TIMES     OF    EVENTS    OP   TIMES     BETWEEN    EVENTS    REQUESTED? 
C 
C 

38  IF    {II.EQ.04    RETURN 
C 

C         CALCULATE    T  IMES     BETWEEN     EVENTS 
C 

S    =    TIMES  (1  ) 

TiMES(l)     =   TI  VESd  )-EL 

N5    =    N(5) 
C 
C 

CO    39    1=2, N5 

SI     =    TIMES( I) 

TIMES(  I)     =    TiyES(  I  )-S 

S    =    SI 

39  CONTINUE 
C 

RETURN 

END 
C         SUBROUTINE    NHPP 2    SIMULATES    A    NON-HO^CGENEOUS 
C         PCISSCN    PROCESS    WITH    A    LOG-LINEAR    INTENSITY 
C         (RATE)    FUNCTION 


C 


C 


SUBROUTINE    NHFP2     (IS , EL , ER, A,  Al ,  II ,N  ,  lER ) 
DIMENSION    T(5C00) 
COMMON    /HOLD/    T 


CALL    OVFLOW 
C 

C         INITIALIZE    VARIABLES 
C 

lER    =    0 

TINT    =    ER-EL 

A    =    EXP(A+A1*EL) 
C 
C       IS    THE    POISSON    PROCESS    HQMGGENEGUS? 


IF     (Al.EQ.O.)     GO    TQ    3 

PAR    =    (A*{EXP(TINT*A1)-1.) ) /Al 

IF     (Al.GT  .0.  )    GO    TO    1 

IFLAG    =    3 

GG    TO    2 

/4    =     A«EXP  {TINT*A1  ) 

Al     =    -Al 
iFLAG    =    2 
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C         COMPUTE   PA.RAMETERS    OF    BOTH    POISSON    RANDOM    VARIABLES 
C 

2  THETA    =  -A/Al 
GO    TO    4 

C         CO'^PUTE    RATE    AND    SCALING    PARAMETERS     FOR    HOMOGENEOUS 

C  PCISSON    PROCESS 

C 

3  PAR    =    TINT*A 
I^LAG   =    1 
41NVRS    =    1  ./A 

C 

C    COMPUTE  NUM3ER  OF  EXPONENTIAL  VARIATES  REQUIRED 


4    NMAX    =    PAR+6.*SQRT(PAR) 
C 

C  IS    This    A    HGMGGENEOUS    POISSON    PROCESS? 

IF  (IFLAG .EQ.l)  GO  TO  17 
C 

C    GENERATE  REALIZATION  ON  POISSON  (THETA)  VARIATE 


5    CONTINUE 

Ci^LL    PVAR     {  IS, THETA, M) 
If     (M.EQ^O)     GC    TO    7 

C 

C         CALCULATE    TIMES    OF    EVENTS 


CALL    SEXPON     (IS,T,NMAX) 

B    =    -Al 

V    =    0. 

JMA>    =    NMAX+1 


C 


c 
c 

DO    6    I=1,JMAX 
C 

C  HAVE    NUMBER    OF    EVENTS    EXCEEDED    THE    VAXIMUM    NUMBER 

C         TFAT    THE    ARRAY    CAN    HOLD? 
C 

IF     (I.GT.NMAX)     GO    TO    3 

V    =    V+T  (I  )/(  (M-I-»-l  )*P) 

IF     (V.GT.TINT)    GO    TO   9 

T( I)    =    V 

IF     (i.EQ.M)     GO   TO    10 

6  CONTINUE 
C 

C 

C         NO    EVENTS    OCCURRED 

C 

7  N    =    0 
RETURN 

C 

C  TCO    MANY    EVENTS    FOR    ARRAY •       INCREMENT    ERROR 

C         CODE    AND    TRY    AGAIN 

C 

8  lER    =    IER+1 
GO    TO    5 

C 

C         THE    NUMBER    OF    EVENTS    OBSERVED    TO    OCCUR    IN    THIS 

C         NON-HOMOGENEOilS    POISSON    PROCESS    IS     'N' 

C 

9  N    =    I-l 

IF    (N.EO.O)    RETURN 
GO    TO    11 

10  N    =    M 

11  CONTINUE 
C 

C    IS  THE  RATE  FUNCTION  INCREASING  OR  DECREASING? 
C 

IF  (IFLAG. EG. 3)  GO  TO  13 
C 

C      TIME    REVERSAL    TECFJ^IOUE     IS    NECESSARY 
C         DETERMINE    WHETHER    N    IS    EVEN    OR    ODD 
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ISIG    =    M0D(N,2) 
NLOQP    =   H/2 

Nl    =    N+i 

c 
c 

CJ    12    I=1,NL0CP 
S    =    TCI) 

T(  I)    =    ER-T(N1-I ) 
KNl-I)    =    ER-S 

12  CONTINUE 

C 

IF  (ISIG.EQ.l)    T( NLCaP+l)=ER-T(NLaQP+l) 
C 

C         ARE    TIMES    OE    EVENTS    REQUESTED? 
C 

IF  (li.EQ.O)     FETURN 

GO  TO    15 

13  IF  (II.NE.a)     GO    TO    15 
IF  (£L,EQ.O.l    RETURN 

C 
C 

00    14   I=lfN 

T(  I)    =    EL+T(I  ) 

14  CONTINUE 
C 

RETURN 
C 
C         CALCULATE    TIMES    BETWEEN    EVENTS 


C 

c 
c 


c 
c 


15    S    =    T(l) 


DO     16    I  =2,N 

SI      =    T{I) 

T(I)    =   T( IJ-S 

S   =    S"* 

16  CONTINUE 
C 

RETURN 
C 

C         THE    PCISSGN    PROCESS    IS    HOMOGENEOUS 
C 

17  I    =    1 
U    =    0. 

CALL    SEXPCN     (ISfTjNf'AX) 

18  M    =    U+T(I ) 

IF    (U.GT.PAR)    GO    TO    20 

I    =    I+i 

IF     (I.GT.NMAX)     GO    TO    19 

GO    TC    18 
C 
C         INCREMENT    ERROR    CODE 


19  lER  =  IER+1 
C 

C    TRY  AGAIN  WITH  NEW  STRING  OF  VARIATES 

GO    TO    17 

20  N    =    I-i 
IF    (N  .EQ.O)    RETURN 
IF     (lI.EO.l  )    GO  TO    22 


DO    21    1=1  ,N 
EL    =    EL+AINVRS*T(I ) 
T  (  I  )    =    EL 
21    CONTINUE 


RETURN 
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c 
c 


c 


c 
c 


c 


22  CO    23    1=1, N 

T( I  )    =   T( I)*AINVRS 

23  CONTINUE 


RETURN 

END 
C         SUBROUTINE    PVAR    GENE'.ATES    A    PQISS3N    (  THETA ) 
C  VARIATE,    M,    USING    THE    GAMMA    METHOD 


SUEROUTINE    PV/^R     (IS, THETA, M) 

DIMENSION    T(5000) 

COMMON    /HOLD/    T 

K    =    0 

C    =    16.0 

D    =     .875 

IF     (THETA. LT.C)    GO    TO    2 

GO    TO    5 

U    =    1. 

CTN    =    EXP (-THETA) 

MMA>    =    THETA+6.*S0RT(THETA) 

I    =    1 

CALL    SRAND    <IS,T,MMAX} 

U    =    U^H  I  ) 

IF     (U.LT.CTN)     GO   TO    8 

I    =    I  +  l 

K    =    K  +  1 

IF     (I  .GT.MMAX)     GO    TO    3 

GO    TO    4 

NP    =    INT(0*THETA) 

AN    =    FLCAT(NP) 

CALL    GAMA     (AN,IS,G,1) 

IF     (G.GT. THETA)    GO    TO    6 

K    =    K+NP 

THETA    =    THETA-G 

GO    TG    1 

L    =    THETA /G 

NP    =    NP-i 

CALL    SRAND    (;IS,T,NP) 


DO    7    I=1,NP 
l^     (T(I  ).LT.U)    K    =    K-t-1 
7    CONTINUE 
C 
C 
C 

C  THE    VALUE    M    IS    ASSUMED    BY    THE    POISSCN     (THETA)    V/5RIATE 

C 

3     1^    =    K 
RETURN 
END 
C  SL8RQLTINE    REJECT    GENERATES    AN    ORCEREC    SAMPLE 

C         OF    GIVEN    SIZE    FROM    THE    SECOND    COMPONENT 
C  OF    THE    ORIGINAL     INTENSITY    FUNCTION 

C         USING    A    REJECT! CN-ACCEPTANCE    ALGORITHM 


SUePCUTINE    REJECT     ( IS , EL, CST AR, PV EC, G, L ) 

DIMENSION    Vi5C0),     PVEC(5) 

DIMENSION    T(5C00) 

COf'VGN   /HCLD/    T 

L20    =    L*10 

IF     (L20.GT.500)    L20=500 

Li    =    L+1 

K    =    1 

J    =    0 

CALL    SRAND    (IS,V,L20) 


00    2    1=1,  L2  0 
J    =    J+1 
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c 
c 


c 
c 

c 


c 
c 


c 
c 


c 
c 
c 


T(K)    =   Q*V(  J)+EL 

J    =    J+i 

if     (V(J).LT  .CALC(PVEC,T(K))/CSTAR )    K=K  +  1 

IF     (K.EQ. LI)     GG    TO    3 

IF    (J  .GE.L20-I)    GO    TO    1 

CONTINUE 


IF     (K.LT.  L)    GC    TO    1 
3    CALL    PXS3RT    ( T,l,L) 

RETURN 

END 
SUBROUTINE    COLATE    SUPERPOSES    TWO    OROEREO 
EVE.MT    STREANS     OVER    THE    SAME    INTERVAL 

SUBROUTINE    COLATE     {  LBGN  ,LFI  N  ,  JBGN  ) 

Oi^ENSION    TiMES(5000),     T(5000) 

COMMON  /MIKE/  TIMES/HOLD/T 

I  =  JBGN 

J  =  I 

K  =  LBGN 

1  IF     (TIMESd  ),LT.T1MES(K  ))    GO    TO    2 
T(  J)     =    TIMES(  K) 

J    =    J  +  1 

K    =    K+1 

IF     (K.GT.LFIN)    GO    TO    3 

GO    TO    I 

2  T(  J)    =    TIMES  (I) 
J    =    J+1 

I    =    I+l 

IF     (I  .EQ.LBGN)    GO    TO    5 

iO    TO    1 


3    JI    =    LBGN-1 


CO    4    N=I, II 

T(  J)     =    TI  MES(N) 

J    =    J+1 

4  CONTINUE 

RETURN 

5  CONTINUE 

DO    6    N=K,LFiN 
T(N)    =    TIMESC  i\) 

6  CONTINUE 


RETURN 

END 
SUBROUTINE    HELP    EVALUATES    THE    INTEGRATED    INTENSITY 
FUNCTION    OVER    THE    INTERVAL    (EL,ER) 


C 
C 


SUB 
OOU 

z  = 

Y  = 

AA 
BB 
CC 

cc 

If= 

Ql 
02 
XX 
RET 
CC 
XX 
RET 
END 
FUNCT 
THE  D 


1 


-J.OUTINE    HELD    (  A  »  Al,  A2»  EL  ,ER,  XX  ) 
BLE    PRECISION    MMCAW,e3tAA 

SQRT( ABS( A2) ) 

{A1«Z  )/  (2.*A2) 
=    Z*EL+Y 
=    Z=^ER+Y 

=    A-Ai=!'Al/(4.*A2) 
=    £XP(CC)/Z 
(A2.LT.G.)    GO    TO    1 
=    CEXP  (AA*=*2  )*^^MOAW(  AA  ) 
=    DEX? (BB**2)*MMCAW(BB) 
=    CC*{  Q2-gi) 
LRN 

=    CC*. 8862259 
=    CC*(DERF  (BB)-OERF(AA)  ) 
URN 


ION    CALC    EVALUATES    THE    SECOND 
ECOMPOSED    INTENSITY    FUNCTION 


COMPONENT    OF 
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c 
c 


FOR    ANY    INPUT    VALUE. 

FUNCTION  CALC  (PtABSA) 

DIMENSION    P(5) 

X    =    P(i  )^■P{2)*ABSA  +  P  (3)*ABSA**2 

XX    =    P(4)  +P(  5)*ABSA 

CALC    =    EXP(X)-EXP (XX) 

RETURN 

END 
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